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VARIATIONS  IN  BODY  WAVE  MAGNITUDE 
FOR 

YUCCA  FLAT,  NEVADA  EXPLOSIONS 


JOHN  F.  FERGUSON 


ABSTRACT 


This  paper  is  an  investigation  of  teleseismic  body  wave  amplitude 
anomalies  caused  by  near  source  scattering.  A  profile  in  the  central 
portion  of  Yucca  Flat  was  selected  for  two  dimensional  seismic  response 
calculation.  The  subsurface  control  on  the  model  structure  is  excellent 
and  variation  in  the  strike  direction  is  minimal.  There  are  thirteen 
nuclear  shots  arrayed  on  the  profile  in  a  reasonably  uniform  manner. 

A  network  of  six  seismic  stations  was  selected  and  the  P-wave  amplitude 
was  calculated  at  five  frequency  values  uniformly  spaced  between  0.8 
and  1.2  Hz  for  each  shot-station  pair.  These  spectral  values  were 

A 

reduced  to  magnitude  by  Bached  mj,  estimator.  The  theoretical  magni¬ 
tudes  can  be  compared  to  the  observed  magnitudes  which  have  been 
adjusted  to  a  common  magnitude  by  a  scaling  relationship.  The  resulting 
pattern  of  magnitude  deviation  compares  very  favorably  with  observations 
by  Alewine  (personal  communication).  The  best  test  of  this  prediction 
would  be  to  obtain  seismograms  from  the  shot  and  station  pairs  used 
in  this  report  and  compute  for  each  pair. 


INTRODUCTION 


A  systematic  variation  in  body  wave  magnitude  has  been  discovered 
at  Yucca  Flat,  Nevada,  through  an  analysis  by  Alewine  (personal 
communication)  of  the  yield-magnitude  relationship  for  nuclear 
explosions.  Comparison  of  the  expected  magnitude  for  explosions  of 
known  yield,  derived  from  scaling  laws  similar  to  those  of  Mueller 
and  Murphy  (1971)  and  Murphy  and  Mueller  (1971),  and  observed 
magnitudes  at  teleseismic  distances  reveals  a  decrease  in  magnitude 
from  west  to  east  of  over  0.3  magnitude  units.  The  pattern  of 
deviation  is  correlative  with  the  local  geologic  structure,  which 
is  a  Cenozoic  graben  with  an  axial  trend  roughly  north  to  south. 

The  basin  is  filled  with  low  velocity  alluvium  and  volcanic  rocks, 
which  contrast  with  the  surrounding  higher  velocity  Paleozoic 
sediments.  Previous  studies  by  Ferguson  (1981,  1982)  and  Hart, 

£t  ^1  (1979)  support  the  hypothesis,  that  the  amplitude  deviation 
is  due  to  scattering  in  the  near-source  structure. 

The  previous  results  discussed  by  Ferguson  (1981,  1982)  were 
limited  in  that  both  the  wavefield  and  structure  were  modelled  in 
two  dimensions.  The  two  dimensionality  assumption  is  valid  for  the 
structure,  at  least  in  selected  regions  of  Yucca  Flat,  but  it  limits 
the  location  of  receiver  locations  to  a  great  circle  containing  the 
model  profile.  In  a  two  dimensional  structure' the  scattering  inter¬ 
action  is  limited  to  the  model  plane,  and  good  qualitative  agreement 
with  observations  was  found  by  averaging  over  incidence  angle. 

Attempts  to  model  specific  shots  have  been  thwarted  by  an  inability 
to  recreate  the  station  network  and  averaging  process  used  to  reduce 


Che  magnitude  observations.  Considerable  bias  results  from  azimuthal 
effects  and  receiver  site-dependent  effects.  The  experiment  presented 
in  this  paper  endeavors  to  overcome  both  of  these  obstacles.  A 
network  of  six  stations  was  selected  along  with  a  profile  through  the 
center  of  Yucca  Flat.  The  stations  are  representative  of  those  used 
to  monitor  nuclear  explosions  (Bache,  personal  communication),  and 
the  profile  was  selected  on  the  basis  of  structural  uniformity  and 


shot  distribution. 


THE  MODELING  TECHNIQUE 


I 


The  model  is  defined  to  have  two  dimensional  structural  variation 
with  homogeneous  regions  of  variable  vertical  thickness  separation  by 
non-planar  interfaces.  A  numerical  approach  is  necessary  to  solve 
the  wave  equation  in  models  of  this  type.  All  calculations  performed 
for  this  study  make  use  of  the  Aki-Larner  method  (Aki  and  Larner,  1970) 
as  described  in  Ferguson  (1981).  The  solution  to  the  wave  equation 
is  expanded  in  a  series  of  plane  waves  within  each  homogeneous  region. 
At  each  interface  the  boundary  conditions  of  continuous  normal  dis¬ 
placement  and  stress  are  satisfied  in  at  least  squares  sense  at  a 
finite  set  of  points.  In  this  manner  the  frequency-wavenumber  spectrum 
of  each  wave  type  (P,  SV,  or  SH-waves)  is  produced  in  each  layer.  The 
derivation  of  the  method  is  contained  in  Larner  (1970)  and  Bouchon 
and  Aki  (1977). 

A  large  set  of  equations  must  be  solved  at  each  interface  for 
each  frequency  component.  The  computational  effort  is  substantial 

even  for  modern  computing  machines  and  is  the  primary  limitation  in _ 

the  large  scale  application  of  the  technique.  Due  to  this  fact,  only 
two  dimensional  structures  are  practical  although  the  wave  field  is 
specified  in  three  dimensions.  Another  factor  of  some  importance  is 
the  Rayleigh  ansatz  error  which  results  from  a  incomplete  description 
of  the  wave  field  in  the  vicinity  of  the  boundary.  '  This  can  be 
overcome  by  use  of  the  extended  boundary  condition  of  Waterman 
(Bates,  1980),  however  the  present  program  does  not  contain  this 
improvement.  The  boundary  errors  in  normal  diplacement  and  traction 


are  carefully  monitored  and  have  not  been  observed  to  be  a  problem.  The 
statistics  of  the  error  magnitude  in  displacement  and  traction  were 
computed  and  a  histogram  plotted  for  each  interface  calculation.  In 
addition  a  plot  of  error  magnitude  versus  horizontal  location  was  plotted 
so  that  correlation  With  interface  shape  could  be  detected.  A  discussion 
of  many  practical  considerations  for  basin  response  calculations  can 
be  found  in  Bard  and  Bouchon  (1980a,  1980b). 

The  source  contribution  to  the  wavefield  can  be  applied  in  two 
ways;  the  source  spectrum  can  be  added  to  the  response  directly  or 
the  reciprocity  theorem  can  be  used  to  generate  the  far  field  response 
from  that  of  an  incident  plane  wave  (Bouchon,  1976).  The  second 
approach  was  chosen  for  this  study  due  to  the  large  saving  in  computer 
time  that  was  anticipated.  Essentially  each  station  defines  a  plane 
P-wave,  with  azimuth  and  incidence  angle  appropriate  to  the  ray  path 
from  Yucca  Flat  to  the  station.  The  P-wave  response  at  the  station 
from  any  number  of  shot  locations  within  the  model  can  be  obtained 
from  the  single  response  calculation  for  the  incident  wave. 

The  source  must  be  specified  in  the  frequency  domain.  In  the 
calculations  presented  here,  source  was  represented  by  the  spectral 
seismic  moment  for  a  Heavyside  unit  step.  The  more  realistic  Von  Seggern 
and  Blandford  (1972)  source  function  was  also  tried,  but  the  resulting 
pattern  of  magnitude  deviations  was  unaffected. 


THE  EXPERIMENTAL  PLAN 


The  locations  of  the  six  seismic  observatorys  simulated  in  these 
calculations  are  shown  in  Figure  1.  They  were  chosen  to  be  repre¬ 
sentative  of  those  stations  used  for  routine  nuclear  explosion 
monitoring.'  Figure  2  is  a  map  of  the  Yucca  Flat,  Nevada,  area,  with 
the  basin  outlined  in  bold  lines  and  the  structural  contours  on  top 
of  the  Paleozoic  shown  as  light  lines.  The  structure  is  derived  from 
seismic  reflection,  gravity  and  borehole  data  in  Ferguson  (1981). 

The  selected  profile  is  drawn  from  A  to  A'  and  the  locations  of 
explosions  within  1600m  of  the  line  are  presented  as  dots.  The 
bearings  of  the  six  stations  are  depicted  on  the  map  border.  The 
station  locations  are  given  in  Table  I  and  data  on  the  shots  are 
tabulated  in  Table  II.  With  these  shot-station  pairs  it  is  possible 
to  compare  the  results  of  this  numerical  experiment  directly  to 
observation  without  the  intermediate  averaging  process  associated 
with  the  reduced  magnitude  deviations. 

Due  to  the  expense  involved  in  synthesis  of  time  domain  signals 
when  using  a  frequency  domain  technique,  it  was  decided  to  use  a 
frequency  domain  magnitude  formula  to  compare  the  theoretical  and 
observational  data.  In  addition  it  has  been  demonstrated  that 
spectral  based  magnitude  estimates  are  more  stable  and  reliable 
than  the  traditional  time  domain  methods.  The  Bache  (1981) 
estimator,  m^,  was  chosen  for  this  application.  This  technique 
was  designed  to  use  spectral  amplitudes  produced  by  the  MARS 
program.  MARS  is  based  on  an  analytic  signal  decomposition  and 


A  STATION 
0  NEVADA  TEST  SITE 


Figure  1.  This  map  displays  the  location  of  Yucca  Flat,  Nevada  as 
a  hexagon  in  the  center  and  the  locations  of  the  six  seismic  stations 
as  triangles.  The  lines  of  latitude,  longitude  and  range  are  graduated 
in  ten  degree  units. 


Estimates  from  Dahlman  Israelson  (1977),  except  Bilby  which  was  announced 


produces  a  rather  smooth  spectral  estimate.  Bache  selected  five 
frequency  samples  between  0.8  and  1.2  Hz  and  smoothed  these  by 
fitting  a  second  degree  polynomial  by  least  squares.  The  log¬ 
arithm  of  the  polynomial  evaluated  at  1  Hz  is  used  to  estimate 
the  magnitude.  The  paths  between  each  station  and  all  shots  are 
essentially  the  same,  so  that  further  corrections  for  path  effects 
are  unnecessary  because  only  the  difference  between  shots  or 
stations  is  all  that  is  required. 

The  basin  was  modelled  by  four  layers  representing  dry 
Quaternary  alluvium,  dry  Tertiary  volcanic  rocks,  saturated  volcanic 
rocks  and  Paleozoic  basement,  which  is  mostly  limestone,  quartzite 
and  other  high  velocity  sedimentary  rocks.  The  properties  of 
these  materials  were  selected  from  an  analysis  of  well  logs  and 
are  presented  in  Table  III. 

In  addition  to  the  velocities  a  Q  of  50  was  specified  for 
the  entire  model.  This  is  necessary  to  insure  numerical  stability 
and  is  also  within  the  expected  Q  range  for  this  region.  The  depths 
of  each  of  the  three  interfaces,  as  determined  from  borehole  and 
the  gravity  model  of  Ferguson  (1981),  and  shown  in  Figure  3.  The 
gravity  model  was  produced  by  a  three  dimensional  inversion 
technique  and  is  indicated  by  a  dashed  line.  The  solid  line  is 
the  interface  actually  used  in  the  calculations.  The  interface 
depths  from  well  logs  are  projected  orthogonally  into  the  profile 
and  are  displayed  on  the  cross  section.  It  is  seen  that  the 
structure  in  this  part  of  Yucca  Flat  is  satisfactorily  controlled 
by  subsurface  information.  Reference  to  structure  contours  on  the 
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top  of  the  Paleozoic  in  Figure  2  provides  support  for  the  assumption 
that  the  basin  structure  is  uniform  along  strike  for  several  kilo¬ 
meters  on  either  side  of  A-A'  (approximately  the  first  Fresnel 
zone ) . 

The  wave  field  is  expanded  in  35  discrete  horizontal  wavenumber 
components  in  the  direction  normal  to  strike.  The  response  for 
incident  P-waves  from  each  of  the  six  stations  were  computed  at 
frequencies  of  0.8  to  1.2  Hz  at  0.1  Hz  intervals.  The  shot  locations 
are  taken  to  be  those  of  the  13  shots  in  Table  II  projected  along 
strike  into  the  profile.  The  amplitude  responses  were  reduced  in 

A 

the  manner  of  the  Bache  m^  estimator  for  a  total  of  78  magnitude 


estimates. 


RESULTS 


The  statistics  for  the  magnitudes  at  each  shot  point  are  printed 
in  Table  IV.  The  distribution  has  a  nearly  coihcident  mean,  median 
and  midrange  for  each  shot  and  is  thus  nearly  symmetric.  The  range 
is  0.3  to  0.4  mt,  units  and  the  standard  deviation  is  near  0.1 
units.  There  is  an  undetermined  additive  constant  for  each  shot, 
which  is  proportional  to  the  moment  and  hence  the  yield  of  the 
explosion.  The  mean  of  each  shot  and  limits  corresponding  to  one 
standard  deviation  are  plotted  as  a  function  of  horizontal  position 
in  Figure  3.  Also  plotted  in  the  same  figure  are  the  actual  observed 
m^  deviations  (Alewine,  personal  communication).  The  magnitude 
level  was  adjusted  so  that  the  theoretical  and  observed  values 
matched  on  the  average  over  the  entire  data  set.  The  pattern  of 
west  to  east  decrease  in  magnitude  is  fit  extremely  well  by  the 
calculations  in  spice  of  the  in  station  network  and  magnitude 
formula  used  in  the  calculation  were  different  from  that  used 
to  find  the  observed  magnitudes  (Alewine,  personal  communication). 

It  is  interesting  to  note  that  the  dispersion  of  the  estimates 
increases  to  the  east  as  the  deviation  decreases.  With  the  in¬ 
creasing  dispersion  a  large  bias  due  to  network  differences  would 
be  expected  and  may  be  observed  in  the  fit.  The  only  shot  to  fall 
outside  of  the  one  standard  deviation  limit  is  Bilby.  This  is  the 
earliest  and  largest  shot  and  may  be  effected  by  a  network  bias 
or  diffidences  in  the  yield  scaling  laws. 

Two  other  calculations  were  performed  with  the  same  parameters 
and  model,  but  different  source  distributions.  The  sources  were 


4.  Table  IV 

Deviation  Statistics  of  for  Each  Shot 

Name _ Mean _ Median _  Midrange _ Standard  Deviation _ Range 


Sandreef 

0.076 

0.063 

0.057 

0.133 

0.365 

Rummy 

0.045 

0.004 

0.005 

0.145 

0.343 

Algodones 

0.091 

0.105 

0.054 

0.102 

0.252 

Wagtail 

0.185 

0.170 

0.193 

0.080 

0.203 

Tan 

0.178 

0.126 

0.210 

0.081 

0.209 

Oscuro 

0.225 

0.209 

0.231 

0.103 

0.306 

Keelson 

0.233 

0.201 

0.262 

0.096 

0.289 

Escabosa 

0.231 

0.194 

0.266 

0.093 

0.275 

Piranha 

0.178 

0.154 

0.211 

0.091 

0.264 

Bilby 

0.145 

0. 10  L 

0. 191 

0.090 

0.252 

Shaper 

0.038 

0.035 

0.040 

0.264 

0.766 

Tijeras 

-0.055 

-0.159 

0.061 

0.189 

0.496 

Monero 

-0.140 

-0.190 

-0.059 

0.189 

0.523 

The  levels  were  adjusted  by  an  arbitrary  constant  to  fit  the  observations,  which  have 
been  reduced  to  a  common  yield.  The  statistics  are  based  on  the  distributions  for  the 
six  stations  in  figure  1.  The  shots  are  ordered  according  to  occurrence  in  figure  3 
from  west  to  east. 
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Figure  3.  The  lower  part  ot  this  figure  shows  the  geologic  structure 
on  profile  A  to  A'  in  figure  2.  The  triangles  are  depths  to  the 
Tertiary  volcanics  and  the  circles  are  depths  to  the  Paleozoic  sediments, 
both  as  determined  in  boreholes.  The  open  circles  are  the  projections 
of  the  shot  locations  into  the  profile.  The  upper  part  of  the  figure 
illustrates  the  correlation  of  magnitude  deviation  with  the  structure. 

The  black  circles  were  calculated  for  this  report  for  the  shot  locations 
shown.  The  dashed  lines  are  +  la  as  determined  for  the  network  of 
figure  1.  The  triangles  are  actual  observations  from  Alewine  (personal 
communication) . 


arrayed  on  planes  jusc  below  the  water  cable  at  550m  and  600m  depths. 
There  was  no  significant  effect  on  the  pattern  in  either  case  when 
compared  to  the  source  distribution  in  Figure  3. 

The  convergence  of  the  Fourier  series  in  wavenumber  was  tested 
by  a  run  with  53  rather  than  35  wavenumber  samples  for  station  6. 

The  results  were  very  similar  to  the  original  calculation  so  that 
it  is  concluded  chat  satisfactory  convergence  is  obtained  with 
the  shorter  series.  There  is  a  tendency  to  ill  conditioning  in 
the  boundary  value  equations  with  the  addition  of  the  higher 
wavenumbers  and  the  inhomogeneous  waves  they  represent. 

A  final  calculation  on  a  model  lacking  the  water  table  inter¬ 
face  was  performed.  This  interface  has  the  lowest  impedance  contrast 
of  the  three  originally  present.  The  pattern  of  deviation  completely 
changed  on  the  eastern  end  of  the  profile.  An  increase  from  west 
to  east  in  magnitude  was  found  with  no  appreciable  decline  on  the 
eastern  end  as  observed  in  the  actual  measurements.  Based  on  this 
experience,  it  is  felt  that  the  four  layer  model  is  close  to  the 
minimum  complexity  necessary  to  explain  the  observations. 


CONCLUSIONS 


The  Yucca  Flat  geophysical  model  utilized  in  this  calculation 
was  derived  from  abundant  geologic  and  geophysical  data  and  is  quite 
detailed.  The  Aki-Larner  technique  has  been  effective  in  computing 
the  plane  wave  response  of  two  dimensional  cross  sections  taken  from 
the  geophysical  model.  The  plane  wave  response  was  used  in  a 
reciprocity  relationship  in  order  to  obtain  the  far  field  P-wave 
response  at  a  selected  set  of  six  stations  due  to  a  set  of  thirteen 
explosions  in  the  cross  section.  A  frequency  domain  body  wave 
magnitude  formula  due  to  Bache,  was  used  to  reduce  these  calculations 
to  a  profile  of  magnitude  deviation  versus  horizontal  location  on 
the  model.  These  predictions  agree  with  the  observations  of 
magnitude  deviation  for  the  same  shot  points.  The  best  test  of  the 
predictions  would  be  to  use  the  Bache  magnitude  formula  on  the  actual 
seismograms  from  the  same  shots  and  stations. 
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PART  I 

INTRODUCTION 

In  recent  years  it  has  become  evident  that  lateral  variations  in 
geologic  structure  at  Yucca  Flat,  Nevada,  exert  controls  over  the  far 
field  waveforms  from  explosions  detonated  there.  A  marked  correlation 
was  observed  between  structural  contour,  Bouguer  gravity,  and 
normalized  body  wave  magnitude  maps.  This  correlation  suggested  that 
a  deterministic  geophysical  model  based  on  an  accurately  defined 
geologic  structure  could  be  used  to  predict  far-field  body  wave 
amplitude  variations.  When  the  examination  of  existing  structural 
interpretations  proved  to  be  deficient,  it  was  decided  to  obtain 
modern  seismic  reflection  profiles.  In  addition,  a  reevaluation  of 
the  United  States  Geological  Survey  (USGS)  gravity  data  used  in 
earlier  structural  interpretations  seemed  to  be  appropriate.  Various 
borehole  and  other  geologic  data  were  obtained  from  Los  Alamos 
National  Laboratory  (LANL)  and  published  sources  to  aid  this  research. 

Yucca  Flat,  Nevada,  is  a  graben-like  structure  of  Cenozoic  age. 
The  underlying  strata  are  mainly  sediments  of  Paleozoic  (Pz)  age. 

Major  volcanic  centers  to  the  west  have  introduced  a  thick  sequence  of 
tuff  units  which  are  generally  less  than  20  million  years  old  (Tv). 
Within  Yucca  Flat  and  adjacent  valleys  alluvium  of  Quaternary  age  (Qal) 
covers  the  volcanics  to  a  thickness  of  over  300  m.  Major  normal  fault 
sets  strike  east-west,  northwest-southeast  and  north-south.  The  last 
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group  are  apparently  the  youngest  and  have  a  substantial  strike  slip 
component  of  motion.  This  brief  geologic  description  will  be 
expanded  after  introduction  of  the  geophysical  model. 

The  interpretative  procedure  used  from  the  outset  of  this 
study  required  the  examination  of  as  many  different  types  of  data  as 
possible  in  order  that  the  resulting  model  be  a  complete  explanation 
of  the  basin  phenomena,  consistent  with  all  observations.  In 
particular,  gravity  observations,  borehole  data  and  seismic  reflection 
profiles  were  simultaneously  used  in  the  interpretative  process.  It 
was  also  found  that  the  structure  must  be  modeled  in  three  dimensions 
to  achieve  reliable  results. 

This  report  presents  the  resulting  three  dimensional  geophysical 
model  for  Yucca  Flat.  Results  of  numerical  wave  propagation  experi¬ 
ments  will  be  presented  in  a  later  report  directed  toward  the 
seismological  aspects  of  this  study. 


GEOLOGIC  INFORMATION 


Geologic  maps  produced  by  the  USGS  and  compiled  at  a  1:24000 
(7.5  minute  quandrangle)  scale  were  available  for  the  entire  study 
area.  These  are  map  sheets  GQ  213,  214,  215,  363,  384,  577,  582, 

746,  and  1327.  The  geologic  base  map  upon  which  the  geophysical  data 
are  compiled  is  derived  from  these  maps.  Additional  geologic  maps 
and  data  have  become  available  from  LANL  and  USGS  sources  concerned 
with  the  Nevada  Test  Site  (NTS).  The  Geological  Society  of  America 
Memoir  110,  Nevada  Test  Site,  edited  by  Eckel  (1968)  contains  a 
great  deal  of  pertinent  information. 

Borehole  stratagraphic  and  geophysical  data  have  been  of  par¬ 
ticular  importance.  Most  of  the  borehole  data  were  obtained  from 
LANL;  consequently,  only  limited  data  from  areas  managed  by  Lawrence 
Livermore  Laboratory  (LLL)  or  from  areas  that  are  of  little  interest 
to  either  group.  The  borehole  geophysical  data  consist  of  summaries 
compiled  by  LANL  of  average  properties  for  various  stratagraphic 
units.  In  some  cases  the  original  logs  were  available.  There  are  36 
such  boreholes  summarized,  including  7  which  penetrate  to  the 
Paleozoic,  17  which  penetrate  into  the  Tertiary,  and  12  which  bottom 
in  the  Quaternary  alluvium.  In  addition  we  have  examined  borehole 
gravity  data  for  41  holes  in  Yucca  Flat.  Control  on  the  Paleozoic 
surface  is  provided  by  the  90  holes  shown  in  figure  1.  The  alluvium- 
volcanic  contact  is  controlled  by  a  total  of  166  holes.  These  holes 
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are  not  the  complete  sample  of  borehole  data  at  Yucca  Flat.  The 
numbering  of  these  boreholes  corresponds  to  various  numbered  areas 
defined  on  Yucca  Flat  for  administrative  purposes.  Figure  1  pro¬ 
vides  a  rough  outline  of  these  areas  through  the  grouping  of  the 
holes  themselves.  At  several  places  in  this  report  it  will  be 
convenient  to  refer  to  the  groupings  by  area  number. 


-n . 


Fig.  1.  Location  map  of  boreholes  contacting  the  Paleozoic  Base¬ 
ment  used  in  this  study.  The  numerical  portion  of  the  borehole  de¬ 
signation  corresponds  to  area  designations  used  in  some  sections  of 
this  dissertation.  The  geologic  base  map  will  be  used  in  all  sub¬ 
sequent  maps.  The  dark  shade  represents  Tertiary  volcanic  rocks, 
the  light  shade  represents  Quaternary  alluvium  and  the  white  areas 
are  Paleozoic  outcrop.  The  Climax  stock  in  the  northern  part  of 
the  map  is  stippled.  Nevada  Central  Zone  coordinates  in  units  of 
feet  are  provided. 
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GRAVITY  DATA 

Some  8459  gravity  stations  are  located  within  an  area  22.5  minute 
on  a  side  around  Yucca  Flat.  The  data  are  from  USGS  open  files  and  were 
collected  over  a  twenty  year  period  by  various  people  using  different 
instruments.  Most  of  these  data,  6383  stations,  are  within  Yucca  Flat. 
Of  the  remainder,  351  are  sited  on  Paleozoic  outcrops  and  1725  are  in 
neighboring  basins  and  on  the  volcanic  outcrops  which  border  the  basins. 
The  magnetic  tape  file  containing  this  compilation  specified  location  to 
0.01  minute,  elevation  to  0.1  foot  and  terrain  corrected  Bouguer  anomaly 
to  0.01  milligal.  The  precision  and  quality  of  these  data  are  variable. 
An  accuracy  of  estimate  of  +0.5  milligal  is  assumed  for  the  overall  data 
set,  although  the  newer  data  is  probably  good  to  +0.15  milligal  (Felch, 
1979).  Figure  2  is  a  contour  map  of  this  Bouguer  gravity  data. 


Fig. 

area. 


2.  Complete  Bouguer  anomaly  map  of  the  Yucca  Flat,  Nevada 
The  contour  interval  is  1  milligal. 
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ANALYSIS  OF  THE  GRAVITY  DATA 


The  abundant  gravity  observations  (over  8000  stations)  in  a  prop¬ 
erly  devised  gravity  inversion,  constrained  by  the  more  sparse  geologic 
and  geophysical  data,  should  produce  a  relatively  detailed  geophysical 
model.  The  model  geometry  most  appropriate  in  this  case  is  that  of  a 
basin  of  low  density  fill  overlying  a  higher  density  half  space.  The 
shape  of  the  model's  lower  perimeter  will  be  sought  as  a  function  of 
position  using  the  observed  gravitational  acceleration.  The  upper 
boundary  is  the  earth's  surface,  which  is  known  in  advance.  If  the  den¬ 
sity  distribution  is  also  known,  the  interface  shape  is  uniquely  deter¬ 
mined  (Smith,  1961),  for  an  infinite  number  of  accurate  observations. 

A  finite  and  inaccurate  data  set  will  in  general  not  permit  a  unique 
solution,  but  some  average  model  can  always  be  found.  The  Inverse  gra¬ 
vity  problem  for  model  shape  is  nonlinear  and  can  be  performed  iterati¬ 
vely  by  techniques  such  as  Gauss-Newton  or  successive  approximations. 

There  exist  many-  published  inversion  methods  which  have  differing 
strengths  and  weaknesses  in  a  given  context.  There  is  no  perfect 
inversion  theory  of  universal  applicability.  The  following  criteria 
were  considered  in  the  choice  of  an  inversion  routine: 

1)  Ability  to  recover  a  best  fitting  solution  (according  to  some 
norm) . 

2)  Ability  to  assess  the  accuracy  of  the  solution. 

3)  Ability  to  include  a  priori  geologic  and  statistical 
information. 

4)  Speed  and  computer  memory  requirements  necessary  to  invert 
thousands  of  observations. 

5)  Provision  for  a  three  dimensional  model  and  variable  density 
contrast,  if  necessary. 
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Trial  and  error  methods  were  rejected  as  too  crude.  Initially  some 
thought  was  given  to  vertical  prism  parameterization  and  linearized 
inversion  similar  to  the  technique  of  Burkhard  and  Jackson  (1976). 
Preliminary  experiments  indicated  that  although  this  technique 
satisfied  most  of  the  criteria,  the  speed  and  memory  economy  were 
insufficient  to  be  practical  for  the  large  data  set.  In  addition, 
discrete  parameterizations  tend  to  produce  models  which  are  rather 
dependent  on  the  parameterization  (i.e.,  block  size  and  location).  An 
example  of  a  block  parameterized  model  based  on  these  same  data  is  con¬ 
tained  in  Felch  (1979)  who  applied  the  Cordell  and  Henderson  (1968) 
method  to  the  Yucca  Flat  data. 

An  integral  equation  formulation  of  the  problem  could  be  specified 
which  would  overcome  at  least  some  of  these  problems.  The  interface 
would  be  represented  by  a  continuous  function  of  the  spatial 
coordinates.  A  prototype  program  of  this  sort  was  written  for  two 
dimensional  basin  models  which  worked  quite  well.  Results  of  that 
program  were  presented  by  Goforth,  e£  a_l.  (1979).  However,  this 
experience  indicated  that  computational  limitations  would  rule  out 
a  three-dimensional  implementation  of  that  technique.  The  method 
finally  selected  is  similar  in  some  respects  to  the  integral  equation 
formulation.  It  is  the  Parker-Oldenburg-Huestis  method  (Parker,  1972, 
Oldenberg,  1974,  Parker  and  Huestis,  1974),  which  uses  Fourier  expan¬ 
sions  to  achieve  computational  speed  through  application  of  the  fast 
Fourier  transform.  The  details  will  be  discussed  in  a  later  section. 
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In  a  local  study  the  regional  gravity  anomaly  is  really  part  of 
the  model  specification.  It  is  that  part  of  the  field  which  is  not 
accounted  for  by  the  model  in  question  and  usually  is  arbitrarily 
parameterized.  We  choose  to  specify  the  regional  gravity  by  a  low- 
order  Fourier  trend  surface,  which  is  fit  to  those  gravity  stations 
occuring  on  the  Paleozoic  outcrop  (i.e.,  the  half  space  of  the  model). 
Additionally,  the  estimation  of  this  anomaly  is  made  a  part  of  the 
inversion  process  by  removal  of  the  modeled  basin  gravity  effect  at 
the  borehole  Pz  control  points  within  the  basin  and  the  use  of  the 
resulting  anomalies  in  the  least  squares  surface  fitting.  This  concept 
was  applied  to  Yucca  Flat  by  Felch  (1979).  The  notion  has  been 
extended  in  this  application  to  a  procedure  which  iteratively  improves 
the  regional  gravity  estimate.  Starting  with  a  surface  with  no  bore¬ 
hole  control,  the  residual  anomaly  is  inverted.  The  depth  residuals  at 
the  boreholes  are  used  to  correct  the  stripping  of  the  basin  anomaly,  a 
new  regional  anomaly  surface  is  better  defined.  It  is  possible  that 
long  wavelength  variations  in  density  contrast,  if  improperly 
specified,  will  be  confounded  with  the  regional  anomaly  produced  by 
this  process.  A  simplified  flow  chart  of  the  entire  procedure  is  pre¬ 
sented  in  figure  3.  Each  box  is  a  separate  computer  program  which 
operates  on  various  intermediate  data  files.  The  process  may  be  cycled 
as  many  times  as  necessary  to  achieve  stable  error  statistics. 
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Fig.  3.  Flow  chart  depicting  the  major  elements  of  the  gravity 
interpretation  procedure.  The  rectangular  boxes  represent  processes 
performed  by  computer  programs  and  the  elliptical  boxes  represent 
data  files  both  generated  and  operated  on  by  the  programs. 
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THE  PARKER-OLDENBURG-HUESTIS 
POTENTIAL  INVERSION 

Parker  (1972)  devised  an  expansion  of  the  integral  equation  for 
the  gravitational  (or  magnetic)  potential  in  the  form  of  a  series  of 
Fourier  transforms.  In  the  following  we  define  the  z-axis  to  be  posi¬ 
tive  downward.  The  anomalous  potential  is  assumed  to  be  due  to  a  layer 
of  contrasting  density  and  variable  thickness,  which  vanishes  outside 
of  some  finite  area,  D,  and  is  specified  by 
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The  vector  r  is  the  projection  of  r  onto  the  plane  z  .  Parker  then 
—  o 

expands  the  second  exponential  In  a  Taylor  series  and  exhanges  the 
order  of  summation  and  Integration  resulting  in  the  rather  novel  repre¬ 
sentation  of  the  potential  spectrum  as  an  infinite  series  of  Fourier 
transforms, 

00 
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The  potential  may  be  differentiated  with  respect  to  z  to  obtain  the 
gravitational  acceleration.  Equation  4  can  be  rearranged  slightly  into 
a  remarkable  triplet  of  formulas : 
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Now  5a  states  the  forward  problem  of  finding  the  gravitational 
field  given  the  density  and  layer  shape.  Equation  5b  states  the  linear 


inverse  problem  of  the  determination  of  density  given  gravity  obser¬ 
vations  and  layer  geometry,  and  5c  states  the  nonlinear  inverse  problem 
for  layer  shape  given  the  gravity  data,  density,  and  upper  surface. 

The  geophysical  model  of  a  basin  would  have  fixed  upper  and  variable 
lower  surfaces.  Variation  of  the  upper  surface,  with  a  fixed  lower 
surface,  could  model  an  intrusion.  The  inverse  problems  may  be  solved 
by  successive  approximations  given  some  initial  guess  for  g(r)  or  p(r_). 
In  this  approach  the  correction  to  the  function  guess  is  a  nonlinear 
step  at  each  iteration  and  hence  has  some  desirable  properties  relative 
to  the  linearized  inversion  schemes  currently  popular.  In  addition, 
for  N  observations  the  Fourier  integrals  nay  be  computed  in  N  log(N) 

operations  for  gridded  data,  whereas  the  integral  equation  methods 
3 

require  order  N  operations  for  the  decomposition  of  a  system  of  linear 
equations. 

Notice  that  the  leading,  linear  term  of  the  expansion  in  equation 
5  is  merely  the  continuation  operator  in  the  wavenumber  domain.  To 
first  order  the  method  is  equivalent  to  Peters'  method  (Peters,  1949 
and  Gunn,  1975)  for  the  interpretation  of  gravity  and  magnetic  data, 
which  has  been  in  use  for  nearly  fifty  years.  This  now  classic  tech¬ 
nique  applied  a  downward  continuation  operator  in  the  space  domain  to 
obtain  a  surface  density  distribution  at  depth  called  an  equivalent 
stratum.  For  a  constant  density  contrast  the  density. variations  on  the 
equivalent  stratum  translate  into  topography  on  a  surface  about  the 
depth  of  continuation.  The  linearized  method  requires  that  the 
topography  be  very  small  relative  to  the  depth.  Parker's  technique 
adds  in  the  higher  order  terms  as  necessary  to  provide  convergence  for 


arbitrary  topography. 
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Some  important  results  on  series  convergence  and  other  properties 
are  given  by  Parker  (1972),  Parker  and  Huestis  (1974),  and  Oldenburg 
(1974).  A  discussion  of  this  theory  from  the  standpoint  of  integral 
equations  may  be  found  in  Dorman  and  Lewis  (1974)  and  in  a  similar 
derivation  from  the  Russian  literature  in  Tsirul'skiy  (1968)  and 
Tsirul'skiy  and  Ospischeva  (1968). 

A  rather  serious  problem  in  the  two  inversion  formulas  arises 
from  the  downward  continuation  operation.  Downward  continuation  is  an 
ill-posed  problem;  that  is,  any  small  perturbation  of  the  gravity 
field,  such  as  that  due  to  noise  in  the  data,  may  result  in  an 
unbounded  perturbation  in  the  continued  field.  Seme  form  of  regulari¬ 
zation  or  smoothing  is  essential  to  bound  for  the  solution  of  such 
problems.  This  is  particularly  important  in  nonlinear  problems  when 
the  inverse  mapping  is  approximate.  Parker,  Oldenburg  and  Huestis  have 
used  a  rather  arbitrary  low-pass  filter  to  eliminate  short  wavelengths 
from  the  solution,  but  our  experience  shows  that  this  is  not  a  satis¬ 
factory  approach.  A  better  filter ~derlve¥~Trom~ETie~work  of  TIkonov 
(1963).  A  regularization  is  defined  as 

4)  min  (m)J  =  mi  ft  u<  u  (r)  -  (r,m)y  dr  *■  <*.  |/m  ( v )  l|  J 

The  operator.  A,  defines  the  forward  problem,  where  U(r)  is  the 
observed  data,  and  m(r)  is  the  approximate  function  sought.  The  first 
term  is  the  data  misfit  and  the  second  is  some  measure  of  model 


smoothness.  The  parameter,  a,  will  be  called  the  regularization 


« 


19 

parameter;  It  controls  the  smoothing  in  the  inversion  for  the  approxi¬ 
mate  model.  The  first  error  term  will  never  be  zero  but  can  be  made 
quite  small  with  a  proper  choice  of  a  regularized  inverse  operator. 

In  Tikhonov  eit  al.  (1968)  several  approximate  inverses  for  the  down¬ 
ward  continuation  problem  are  given.  The  function, 

7)  A(k)  *  (r0)]  exp  (I  l!<r  exp  (-Ikl20)); 

was  selected  as  most  suitable  for  this  purpose.  It  has  the  form  of  a 
smooth  low  pass  filter  which  cuts  off  a  faster  rate  than  the  downward 
continuation  operator  grows.  The  regularization  parameter,  a,  is 
choosen  by  forming  a  family  of  approximate  inverses  by  varying  a  in  a 
systematic  way,  testing  each  one  by  forward  calculation  and  computing 
the  rms  error  between  the  theoretical  and  the  observed  accelerations. 

The  simplex  procedure  of  fielder  and  Mead  (1965),  as  presented  by  O'Neil 
(1971),  provides  a  convenient  method  for  this  minimization.  In  one 
dimension  (the  a  dimension)  the  simplex  is  a  line  segment  which 
attempts  to  bracket  the  minimum.  This  line  segment  will  flip  back  and 
forth  and  expand  or  contract  as  necessary  to  obtain  the  minimum. 
Convergence  is  reached  when  the  segment  shrinks  below  a  certain  length. 

Unlike  the  methods  based  on  linearized  solutions  to  the  integral 
equations,  the  Fourier  transform  method  does  not  produce  the  Frechet 
derivatives  as  a  by-product.  This  complicates  the  calculation  of 
Backus-Gilbert  linear  resolving  kernels.  Parker  and  Huestis  (1974) 
offered  an  approximate  resolving  measure  for  the  linear  inverse  problem 
(5b),  but  it  performed  poorly.  The  resolving  widths  were  larger  than 
the  model  even  for  the  highest  variance  solutions.  The  space  domain 
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width  of  the  regularizing  filter  provides  some  information  on  resolu¬ 
tion,  but  it  lacks  the  space  variable  Quality  of  the  Backus-Gilbert 
resolution  measures.  Because  of  the  lack  of  a  linear  mapping  it  is 
also  difficult  to  compute  model  variance.  It  is  conjectured  that  the 
magnitude  of  the  variance  could  be  inferred  by  consideration  of  the 
terms  which  are  suppressed  by  regularization.  Acknowledging  these  dif¬ 


ficulties  the  error  analysis  utilized  here  will  be  based  on  the  errors 
between  the  model  and  known  (borehole)  depths. 


THE  DENSITY  MODEL 


The  geophysical  model  proposed  here  consists  of  a  Paleozoic  half 
space  with  a  lower  density  basin  with  Cenozoic  age  fill  on  top.  The 
density  contrast  could  have  any  manner  of  spatial  variability.  The 
simplest  density  distribution  one  can  assume  in  order  to  model  the 
basin  would  be  a  single  constant  density  contrast  across  the  interface 
between  the  Cenozoic  and  Paleozoic  strata.  However,  in  the  Great  Basin 
it  is  common  to  encounter  density  functions  which  increase  more  or  less 
exponentially  with  depth.  Cordell  (1979)  has  proposed  that  lateral 
variations,  with  density  decreasing  toward  the  basin  center  may  be 
important.  Some  basis  for  the  choice  of  a  density  model  must  be  found 
in  the  observed  density  data. 

For  many  years  a  constant  density  contrast,  between  the  Cenozoic 
and  Paleozoic  rocks,  of  -0.7  gm/cc  has  been  the  standard  model  for 
Yucca  Flat  (Hazelwood,  et.al.,  1963,  Healey,  1966  and  1968).  It  has 
also  been  recognized  that  there  are  departures  from  this  average,  but 
the  importance  of  the  departures  has  not  previously  been  determined. 
Density  log  information  from  34  boreholes  was  used  in  a  simple  sta¬ 
tistical  analysis  of  this  question.  This  information  is  from  a  zone 
extending  from  roughly  the  center  of  the  valley  to  the  eastern  margin 
and  covering  about  the  middle  fourth.  It  is  questionable  how  represen¬ 
tative  of  the  remainder  of  the  valley  these  densities  nay  be.  Table  I 
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summarizes  the  statistics  on  the  geologic  units  Qal,  Tv,  and  Pz.  A  t 
test  for  the  equality  of  the  means  of  the  Qal  and  Tv  population  was 
performed  with  the  conclusion  that  the  means  were  equal  at  the  90%  con¬ 
fidence  level.  The  mean  and  standard  deviation  of  the  combined  popula¬ 
tion  are  1.79  gm/cc  and  0.18  gm/cc.  Figure  4  is  a  plot  of  the  histo¬ 
grams  of  the  Qal  and  Tv  samples  with  the  normal  distribution  parameters 
(mean  and  standard  deviation)  for  the  combined  population  superimposed. 
The  frequency  distributions  are  not  normal  and  are  skewed  towards  the 
low  density  end.  A  major  deficiency  in  this  analysis  is  the  lack  of 
volumetric  weighting  of  the  density  values.  Samples  from  thin  layers 
are  given  equal  weight  with  samples  from  thick  layers.  The  conclusion 
is  that  the  Qal  is  indistinguishable  from  the  Tv  in  density,  and  hence 
only  one  layer  of  fill  with  no  vertical  density  variation  need  be  con¬ 
sidered  for  the  valley  model.  The  quality  of  the  summary  data  from 
density  logs  probably  does  not  warrant  a  more  detailed  analysis. 

Better  data  for  the  purpose  of  density  modelling  can  be  obtained 
from  the  borehole  gravity  logs.  Borehole  gravity  provides  a  better 
average  density  for  the  area  around  the  hole  than  density  logs  which 
may  be  complicated  by  local  borehole  effects;  however,  lateral 
variations  in  structure  may  bias  the  gravity  results.  The  borehole 
gravity  logs  which  start  at  depths  of  less  than  200  meters  and  extend 
to  depths  greater  than  500  meters  were  selected  from  the  41  logs 
available.  The  derived  density  and  the  average  density  were  plotted  as 
a  function  of  depth.  The  latter  function  is  more  useful  because  it 
rectifies  the  data  and  suppresses  the  effect  of  thin,  anomalously 
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Fig.  4.  Histograms  of  the  occurance  frequency  of  density  values 
for  all  stratigraphic  units  in  the  Quaternary  (a)  and  Tertiary  (b) 
rocks  of  Yucca  Flat.  The  data  were  compiled  from  interval  density 
estimates  from  the  central  latitude  portion  of  the  valley.  The  bar 
scale  at  the  top  of  each  graph  shows  the  mean  and  one  standard  devi¬ 
ation  fiducials  for  the  combined  distribution. 
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dense,  layers.  When  the  average  density  curve  changes  slope,  it  indi¬ 
cates  a  profound  density  change.  A  comparison  of  the  Yucca  Flat 
average  density  profile  and  those  of  several  other  western  localities 
in  the  United  States  are  shown  in  figure  5.  The  comparison  illustrates 
the  unusual  nature  of  the  Yucca  Flat  distribution.  The  Yucca  Flat  cur¬ 
ves  are  rather  constant  with  depth  and  show  little  evidence  of  a  sig¬ 
nificant  exponential  increase  in  density  with  depth;  many  curves  even 
show  a  decrease  in  density  with  depth.  The  limits  drawn  on  the  graph 
are  the  mean  and  standard  deviation  from  the  geophysical  log  distribu¬ 
tion.  Although  lateral  variation  is  present  it  is  generally  less  than 
+0.1  gm/cc  and  no  systematic  variation  is  apparent. 

The  density  of  the  Paleozoic  basement  is  less  well  determined. 
Some  samples  have  been  tested  in  the  laboratory  and  some  geophysical 
data  exist.  The  choice  of  2.5  gm/cc  seems  to  be  justified  by  these 
data.  Lateral  variations  are  known  to  exist,  but  we  have  insufficient 
information  to  characterize  them.  These  observations  support  a 
constant  density  contrast  model  of  -0.7  gm/cc.  With  more  data, 
variations  of  +0.1  gm/cc  could  be  included  in  the  model,  but  this  is 
probably  unnecessary  for  this  study. 
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Fig.  5.  (a)  is  a  plot  of  average  density  as  a  function  of  depth 

for  Yucca  Flat,  Nevada  and  (b)  shows  the  same  data  for  several  other 
representative  localities  in  the  Basin  and  Range  province.  All  density 
estimates  are  from  borehole  gravity  measurements  except  those  in 
Hachita  and  Sulphur  Springs  valleys.  The  fiducials  lines  drawn  verti¬ 
cally  on  the  left  graph  are  the  mean  and  one  standard  deviation, 
derived  from  the  Quaternary  and  Tertiary  density  data  used  ih  Table  I. 
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AVERAGE  OENSITY  VS  DEPTH  FUNCTIONS 
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THE  GRAVITY  MODEL 

As  outlined  previously  a  gravity  interpretation  scheme  employing 
the  Parke r-Oldenburg  technique  formed  the  heart  of  the  structural 
interpretation.  A  basin  model  was  estimated  in  conjunction  with  the 
regional  anomaly.  The  density  function  was  specified  a  priori.  Ninety 
boreholes  were  used  as  intrabasin  control  points.  The  gravity  field 
and  the  interface  were  specified  on  a  rectangular  grid  with  the  y-axis 
north  and  the  x-axis  east.  The  grid  is  specified  by  76  x  nodes,  65  y 
nodes,  1600  feet  x  spacing,  2000  feet  y  spacing  and  origin  620000  E  and 
780000  N  in  the  Nevada  central  zone  coordinate  system  (units  are  feet). 
The  gridding  is  required  for  efficient  computation  of  the  Fourier 
integrals  in  the  inversion  process  using  the  fast  Fourier  transform. 
Since  the  density  was  assumed  to  be  homogeneous  in  the  basin  it  was  not 
specified  on  a  grid  as  it  otherwise  would  have  been.  The  gridding  was 
accomplished  by  use  of  Shepard's  bivariate  interpolation  function 
(Shepard,  1968)  combined  with  a  rapid  search  technique  adapted  for  the 
high  data  density  (6383  stations)  encountered  at  Yucca  Flat.  Some 
smoothing  is  inherent  in  the  gridding,  which  helps  to  suppress  aliasing 
and  high  frequency  noise  in  the  grid. 

Proper  determination  of  the  regional  gravity  anomaly  is  of  par¬ 
ticular  importance  at  Yucca  Flat.  Inspection  of  small  scale  Bouguer 
anomaly  maps  of  the  western  United  States,  such  as  Eaton,  et  al., 
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(1978),  reveal  that  the  valley  sits  astride  a  high  ridge  extending 
north  from  the  California  border.  Just  to  the  northwest,  a  large 
volcanic  caldera  complex  produces  a  large,  circular  low  which  trun¬ 
cates  this  trend.  The  locally-defined  regional  anomaly  should  account 
for  this  gravity  ridge. 

The  regional  anomaly  was  defined  by  a  low  order  Fourier  surface  of 
two  harmonics  (25  terms).  The  fundamental  wavelengths  of  100000  feet 
east  to  west  and  220000  feet  north  to  south  are  roughly  twice  the  basin 
dimensions.  This  surface  was  computed  by  least  squares  fit  to  350  gra¬ 
vity  stations,  which  occur  in  areas  of  Paleozoic  outcrop,  and  gravity 
anomalies  corrected  for  the  basin  effect  at  the  90  borehole  locations 
in  Yucca  Flat.  The  data  from  stations  and  holes  were  weighted  so  that 
both  groups  of  data  would  provide  approximately  equal  effects  on  the 
resultant  surface.  The  regional  anomaly  is  displayed  in  figure  6. 

This  map  shows  the  ridge  and  part  of  the  expected  circular  depression. 
The  anomaly  is  gently  convex  upward.  If  it  were  not,  and  linear  gra¬ 
dients  were  employed  as  is  common  practice,  the  depth  would  be  substan¬ 
tially  underestimated. 

The  residual  anomaly  at  each  iteration  was  multiplied  by  an  areal 
rectangular  truncation  window  of  highly  irregular  shape,  which  was 
somewhat  larger  than  the  basin  itself  but  smaller  than  the  total  grid. 
To  the  extent  that  the  regional  separation  removed  the  effect  of  the 
half  space  of  the  model,  a  satisfactory  taper  was  achieved.  In  the 
southeast  corner  of  the  basin  the  Paleozoic  rocks  do  not  outcrop  and 
hence  a  fundamental  assumption  of  the  model  is  violated.  The  residual 


Fig.  6.  The  regional  gravity  anomaly  used  in  the  calculation  of 
the  final  depth  model.  A  low  order  Fourier  trend  surface  was  fit  to 
gravity  stations  on  the  Paleozoic  outcrop  and  estimated  anomalies  at 
the  borehole  control  points  shown  in  figure  1.  The  contour  interval 
is  1  milligal. 
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does  not  taper  and  is  abruptly  truncated  by  the  window.  The  resulting 
Gibb's  oscillations  do  not  seem  to  be  significant  within  the  basin.  In 
future  calculations  a  better  window  will  be  developed.  The  residual 
which  was  invented  for  the  final  model  is  mapped  in  figure  7. 

In  each  iteration  of  the  interpretation  process  (see  figure  3)  the 
error  statistics  at  the  boreholes  were  monitored.  After  six  iterations 
no  real  improvement  in  the  fit  was  obtained  and  the  procedure  was 
terminated.  The  final  model  is  referred  to  as  FIB  and  is  shown  in 
figure  8.  The  gravity  misfit,  as  shown  in  figure  9,  is  acceptably  low 
over  the  deeper  parts  of  the  basin  and  increases  slightly  in  the 
shallow  areas.  The  error  statistics  are  summarized  in  Table  II; 
histograms  of  the  gravity  data  misfit  and  borehole  misfit  are  provided 
in  figure  10.  The  distributions  are  not  Gaussian  and  have  some  very 
significant  outliers.  The  rank  order  statistics,  the  median  and  inter¬ 
quartile  deviation,  are  much  more  reliable  estimates  of  location  and 
scale  in  this  case.  The  model  thus  satisfies  the  data  within  +0.5  mgal 
with  some  negative  bias.  The  borehole  depth  mismatch  is  biased  slightly 
positive  which  would  agree  with  the  data  error.  The  borehole  data  have 
several  known  mislocations  and  some  ambiguity  in  the  definition  of  the 
Paleozoic  contact  is  introduced  by  a  thick  colluvium  on  the  Paleozoic 
rocks.  The  stated  borehole  depths  do  not  always  agree  with  the  modeled 
depths  due  to  this  ambiguity.  It  should  also  be  remembered  that  the 
model  depths  are  really  spatial  averages  and  not  point  measurements  like 
the  boreholes.  The  basement  depth  deviation  would  be  robustly  estimated 
to  be  +  80  m,  which  represents  about  a  15%  error  on  the  average.  These 
statistics  are  in  good  agreement  with  Felch  (1979). 
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Fig.  7.  The  residual  gravity  anomaly,  which  results  from  the 
subtraction  of  the  regional  anomaly  in  figure  6  from  the  Bouguer 
anomaly  in  figure  2.  The  contour  interval  is  1  milligal. 
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Fig.  8.  The  Tertiary-Paleozoic  contact  depth  model  for  Yucca  Flat, 
Nevada,  designated  FIB.  Note  the  deep  low  areas  on  the  east  and  the 
shallow  lows  on  the  west  separated  by  a  horst.  The  contour  interval 
is  100  meters. 


Fig.  9.  The  error  between  the  residual  Bouguer  anomaly  of  figure 
7  and  the  theoretical  anomaly  computed  from  the  model  in  figure  8. 
The  contours  are  in  units  of  0.5  milligal.  The  error  is  less  than 
1  milligal  except  in  shallow  areas  near  the  basin  margins  and  over 
the  horst. 
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TABLE  II 


ERROR  STATISTICS  FOR  MODEL  FIB  YUCCA  FLAT, 


DEPTH  MISFIT  AT  90  30RZH0LES  IN  KM 
MEAN  =  0.031 

STANDARD  DEVIATION  =  0.121 
ROOT  MEAN  SQUARE  =  0.125 
MEDIAN  =  0.026 

INTERQUARTILE  DEVIATION  =  0.033 
QUARTILE  COEFFICIENT  OF  SKEWNESS  =  -0.03 


GRAVITY  MISFIT  AT  1336  GRID  NODES  IN  MGAL 


MEAN  =  -0.18 

STANDARD  DEVIATION  =  1.10 
ROOT  MEAN  SQUARE  =  1.12 
MEDIAN  =  -0.39 

INTERQUARTILE  DEVIATION  =  0.49 
QUARTILE  COEFFICIENT  OF  SKEWNESS  =  -0.69 


NEVADA 


Fig.  10.  (a)  shows  the  frequency  of  depth  error  values  in  meters 

at  the  borehole  locations  on  figure  1  for  the  model  shown  in  figure 
8.  (b)  is  a  histogram  of  gravity  misfit  errors  at  the  grid  nodes 

within  Yucca  Flat  which  are  contoured  in  figure  9.  The  fiducial 
limits  of  one  standard  deviation  about  the  mean  are  placed  above  both 
histograms. 
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This  model  is  only  a  first  effort  and  has  limitations  on  it  which 
are  rather  artificial  at  this  time.  In  the  Yucca  Flat  area  over  6000 
gravity  stations  are  available  and  only  1300  grid  nodes  occur  in  the 
same  area.  Thus,  the  data  grid  could  be  considerably  refined.  About 
100  borehole  gravity  soundings  are  available  in  addition  to  a  selection 
of  density  logs.  These  data  could  be  used  to  construct  a  laterally  and 
vertically  variable  density  function  which  would  undoubtedly  improve 
the  fit  in  some  areas.  In  addition  some  technical  improvements  to  the 
algorithm  especially  in  the  areas  of  windowing,  regional  surface 
fitting  and  model  regularization  criteria  are  possible.  These  limita¬ 
tions  in  the  data  processing  result  in  a  decreased  resolution  and  in¬ 
creased  smoothness  compared  to  what  is  actually  possible  using  the 
available  gravity  data. 
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SEISMIC  REFLECTION  PROFILES 

Three  seismic  reflection  lines,  one  oriented  north  to  south 
designated  N-3,  and  two  east  to  west,  labeled  E-l  and  E-3,  were 
recorded  in  January  1978.  The  locations  of  the  lines  are  shown  in 
figure  11.  Four  Y-900  vibrators  provided  sweeps  of  8  to  32  Hz  for  18 
sec  on  N-3.  Sweeps  of  6  to  32  Hz  for  12  sec  were  used  on  lines  E-l  and 
E-3.  Low  frequency  (4.5  Hz)  geophones  were  used  with  330  foot  group 
intervals  on  N-3  and  220  foot  group  intervals  on  the  east  to  west 
lines.  All  lines  had  24-fold  coverage.  The  north  to  south  line,  which 
was  oriented  parallel  to  the  strike  in  an  area  with  generally  thin 
alluvial  cover,  was  intended  to  recover  deep  crustal  reflections,  if 
possible.  A  more  detailed  description  of  the  seismic  reflection  profi¬ 
les  is  contained  in  Goforth,  ex  al.  (1979).  The  data  processing  will 
be  discussed  at  some  length  in  the  next  section. 


Fig.  11.  The  location  of  the  seismic  reflection  profiles  E-l, 
E-3  and  N-3.  The  borehole  control  points  are  also  plotted.  The 
small  numbers,  by  the  profiles,  refer  to  vibration  point  numbers. 
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PROCESSING  OF  THE  SEISMIC  REFLECTION  PROFILES 
AND  VELOCITY  ANALYSIS 


The  three  seismic  reflection  profiles  available  for  analysis  were 
VIBROSEIS^)  data  and  were  initially  cross  correlated  with  the  source 
waveform  to  recover  zero  phase  wavelet  seismograms.  The  records  were 
static  corrected  for  elevation  to  a  datum  of  4300  ft.  and  first  arriv¬ 
als  were  muted.  VELAN^D  velocity  scans  were  performed  at  numerous 
intervals  on  all  three  lines.  VELAnQD  analysis  is  a  standard  velocity 
scan  program  which  computes  trial  normal  moveout  corrections  and  tests 
for  coherence  within  certain  time  windows.  High  coherence  indicates 
that  the  rms  velocity  (approximated  by  the  stacking  velocity)  has  been 
correctly  determined  and  thus  the  velocity  as  a  function  of  depth  is 
determined.  There  are  the  critical  assumptions  of  horizontal  layering 
and  lateral  homogeneity  inherent  in  this  procedure.  Departures  from 
these  assumptions  result  in  over-estimates  of  the  rms  velocity.  For 
homogeneous  layers,  dip  alone  will  produce  an  over-estimate  propor¬ 
tional  to  the  cosine  of  the  dip  angle.  Lateral  heterogeneity  within 
the  limits  of  the  geophone  spread  along  with  diffraction  phenomena 
result  in  even  greater  errors.  This  is  complicated  by  the  fact  that 
reliable  velocity  estimates  require  a  large  geophone  spread  in  order  to 
sample  the  large  moveout  on  the  tail  of  the  travel  time  curve.  Within 
Yucca  Flat  strong  lateral  variation  due  to  faulting  resulted  in  very 
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poor  velocity  estimates  at  depth;  the  shallow  velocity  structure,  down 
to  approximately  500  m  was  adequately  determined.  Many  of  the  VELAN(r) 
profiles  were  unusable,  and  those  that  were  used  produced  structural 
relief  estimates  which  were  inconsistant  with  the  gravity  and  borehole 
analysis.  The  seismic  sections  were  migrated  by  a  finite  difference 
scheme  after  stacking  and  are  displayed  in  figures  12,  13,  14. 

The  lateral  velocity  variations  result  in  a  degraded  common  mid¬ 
point  stack  (CMP)  and  subsequent  migration  of  the  stack  is  also  in 
error.  The  solution  to  this  difficulty  lies  in  migration  of  the 
original,  unstacked  data  and  possibly  the  use  of  migration  velocity 
determination  analogous  to  the  normal  moveout  velocity  (the  normal 
moveout  correction  and  stack  is  the  correct  migration  of  data  observed 
over  horizontal  stratification).  Several  approaches  to  this  problem 
have  been  outlined.  Satlegger  (1975)  and  Dohr  and  Stiller  (1975)  pro¬ 
posed  migration  velocity  determination  and  Doherty  and  Claerbout  (1976) 
offer  a  procedure  for  approximate  migration  before  stack.  These  pro¬ 
cesses  are  complicated  by  the  fact  that  migration  calculations  are 
relatively  more  costly  to  perform  than  normal  moveout  corrections  in 
that  the  better  algorithms  are  based  on  wave  theory  rather  than 
geometrical  ray  theory.  Schultz  and  Claerbout  (1978)  point  out  that 
CMP  gathers  and  stacks  are  generally  not  realizable  wave  fields  espe¬ 
cially  where  there  are  significant  lateral  velocity  variations.  They 
propose  to  migrate  stacks  of  constant  ray  parameter  (slant  stacks).  In 
this  process  a  number  of  beams  would  be  formed  and  each  beam  migrated 
separately,  a  very  costly  procedure.  The  E-3  section  was  processed  in 


Fig.  12.  This  section  is  a  migrated,  common  midpoint  stack  of  the 
east-west  seismic  line  E-l.  The  vibration  points  are  numbered  across 
the  top.  The  total  two-way  travel  time  is  4  seconds. 


Fig.  13.  This  section  is  a  migrated,  common  midpoint  stack  of  the 
east-west  seismic  line  E-3.  The  vibration  points  are  numbered  across 
the  top.  The  total  two-way  travel  time  is  4  seconds. 


Fig.  14.  This  section  is  a  migrated,  common  midpoint  stack  of  the 
north-south  seismic  line  N-3.  The  vibration  points  are  numbered  across 
the  top.  The  total  two-way  travel  time  is  4  seconds. 


a  manner  similar  to  the  method  described  by  Doherty  and  Claerbout 
(1976).  Each  normal  moveout  corrected  CMP  gather  was  migrated  indepen¬ 
dently  with  a  rather  arbitrary  choice  of  velocity  function,  a  new  velo¬ 
city  analysis  was  performed,  and  the  migrated  gathers  were  stacked  with 
the  new  velocity.  This  processing  was  more  expensive  (24  times)  than 
conventional  processing,  but  not  as  expensive  as  migration  of  common 
source  gathers  would  have  been.  The  "migration  before  stack"  section 
is  shown  in  figure  15.  The  new  velocity  analysis  improved  the  deeper 
velocity  estimates  by  over  ten  percent. 

The  choice  of  a  velocity  function  suitable  for  use  in  depth  esti¬ 
mation  was  complicated  by  the  lack  of  good  velocity  scans  and  the  spar¬ 
sity  of  deep  velocity  logs.  It  was  not  possible  to  construct  laterally 
variable  velocity  models;  however,  reasonable  results  were  obtained 
with  a  single  velocity  function.  A  comparison  of  interval  velocity 
data  from  areas  3,4,6,  and  7  is  made  in  figure  16.  The  velocity  func¬ 
tion  defined  by  borehole  UE6E  is  close  to  the  median  and  is  the  deepest 
available  hole.  An  ~uphole~ velocity- survey" conducted  in_UE10BD  by  LANL 
and  LLL  agrees  quite  well  with  that  in  UE6E.  Some  holes  in  Area  7  and 
Area  3  have  much  lower  velocities  down  to  a  600  m  depth  and  some  Area  7 
holes  show  rather  higher  velocities  at  shallow  depth.  Figure  17  is  a 
comparison  of  the  UE6E  velocity  function  with  the  orginal  VELAN^D  and 
"migration  before  stack”  velocity  functions.  The  velocity  analyses 
chosen  are  the  best  available.  They  agree  quite  well  in  the  shallow 
regions,  but  diverge  at  depth.  The  conventional  VELAnOD  would  produce 
errors  of  100  m  in  depth.  The  UE6E  velocity  function  was  chosen  as  the 


Fig.  15.  Seismic  reflection  profile  E-3.  This  is  the  same  data 
shown  in  figure  13  with  different  processing.  In  this  case  the 
common  midpoint  gathers  were  migrated  prior  to  stacking  and  a  new 
stacking  velocity  was  derived  which  was  less  effected  by  lateral 
velocity  variation.  The  vibration  points  are  numbered  across  the  top. 
The  total  two-way  travel  time  is  4  seconds. 


Fig.  16.  Average  velocity  as  a  function  of  depth  determined  from 
interval  velocities  at  Yucca  Flat,  Nevada.  The  curve  labeled  UE6E 
is  from  a  deep  borehole  located  in  the  southern  portion  of  the  basin 
(see  figure  11).  The  UE6E  curve  represents  the  median  choice  for 
the  velocity  function. 


Fig.  17.  The  veolcity  functions  derived  from  stacking  velocities 
and  the  UE6E  interval  velocities.  SP150  on  profile  E-3  is  located 
only  5  kilometers  from  UE6E  in  an  area  of  rather  uniform  structure. 
Notice  that  the  migration  before  stack  function  is  in  better  agreement 
with  the  UE6E  function  at  depth  and  that  little  difference  is  present 
among  the  three  estimates  above  700  m  depth. 
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standard  and  has  been  used  in  all  depth  estimates  from  the  reflection 
data.  This  function  may  need  modification  in  Areas  3  and  7  (profile 
E-3  is  in  Area  3),  but  it  serves  quite  well  in  most  of  the  basin. 
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ANALYSIS  OF  THE  SEISMIC  SECTIONS 
AND  GRAVITY  MODEL 

The  three  seismic  sections,  E-l,  E-3  and  N-3,  have  been  analyzed 
in  an  attempt  to  define  the  most  prominent  geologic  features  of  the 
Yucca  Flat  Basin.  The  seismic  sections  have  been  compared  to  the  grav¬ 
ity  model  and  to  the  borehole  data  with  the  UE6E  velocity  model  pro¬ 
viding  depth  to  time  conversion. 

The  FIB  gravity  model  and  geologic  map  shown  in  figure  18 
demonstrate  several  important  features  of  the  basin.  The  deepest  part 
of  the  basin  is  in  the  southeast.  A  saddle  separates  the  deeper  basin 
from  a  bilobate  northern  basin.  The  saddle  trends  northwest  along  the 
strike  of  several  faults  seen  on  the  eastern  and  western  margins  of  the 
valley.  The  western  boundaries  of  these  deeper,  eastern  basins  will  be 

referred  to  as  the  Carpetbag  Fault  system  as  the  fault  segment  so  named 

is  coincident  with  it.  The  better  known  Yucca  Fault  runs  down  the 
center  of  the  eastern  basins.  Both  the  Yucca  and  Carpetbag  faults  have 
right  lateral  strike  slip  component  of  motion  as  well  as  dip  slip.  A 
narrow  horst,  which  rises  to  within  a  hundred  meters  of  the  surface, 
isolates  a  string  of  shallower  smaller  basins  on  the  west  from  the  main 
eastern  basins. 

The  N-3  line  unfortunately  runs  parallel  to  the  western  basin  and 
slightly  to  the  west  of  its  center.  This  causes  the  reflections 

recorded  to  migrate  from  out  of  the  plane  of  the  section,  thus  pro- 


Fig.  18.  The  three  seismic  profiles  superimposed  on  a  map  of  the 
Tertiary-Paleozoic  interface  and  normal  faults  from  the  geologic  maps. 
Notice  the  Yucca  Fault  in  the  middle  of  the  basin  and  also  the  north¬ 
west  trending  faults  in  the  northwest  and  east  central  portions  of  the 
map,  which  align  with  the  saddle  in  mid-basin.  The  contour  interval 
is  100  m. 
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viding  more  limited  resolution  of  vertical  depth.  The  N-3  section  was 
successful  in  recovering  structural  information  from  below  the  basin 
and  possibly  a  mantle  reflection  at  10.7  seconds  two-way  travel  time 
near  vibration  point  140. 

The  two  east  to  west  trending  seismic  lines  provide  very  good 
structural  information.  In  order  to  facilitate  the  interpretation  of 
these  sections  the  FIB  model  was  converted  into  time  sections  along  the 
three  reflection  lines  and  then  overlain  on  the  migrated  common  depth 
point  sections,  as  shown  in  figures  12,  13,  and  14.  Prominent 
reflecting  events  and  faults  delineated  by  diffractions  were  traced 
onto  the  overlays,  as  seen  in  figures  19,  20,  and  21.  In  many  areas 
the  correlation  between  the  gravity  anomaly  and  the  seismic  reflections 
is  quite  good,  but  in  other  areas  the  Paleozoic  contact  does  not  appear 
to  be  a  conspicuous  reflector.  In  these  cases  the  deeper  reflections 
may  be  masked  by  thin,  high  velocity  welded  tuff  units  at  shallower 
depths.  The  horst,  bounded  on  the  east  by  the  Carpetbag  Fault  system, 
is  evident,  as  is  the  Yucca  Fault.  Dip  slip  motion  on  these  faults, 
especially  the  Carpetbag,  accounts  for  the  asymmetric  basin  that  has 
recently  developed  at  Yucca  Flat.  The  dips  of  these  faults  seem  to  be 
somewhat  low  for  normal  faults,  between  45°  and  60°.  Note  that  the 
curved  lines  drawn  on  che  time  sections  (which  map  nonlinearly  into 
depth)  would  appear  as  planes  in  a  depth  section.  West  dipping  faults 
on  the  eastern  end  of  the  lines  have  steeper  dips  on  the  order  of  70° 
and  are  conjugate  to  the  basin  forming  faults.  In  section  E-3  the 


water  table  accounts  for  the  horizontal  reflector  at  0.7  seconds  travel 


Fig.  19.  An  interpretation  of  the  first  second  of  two-way  travel 
time  for  the  E-l  profile  shown  in  figure  12.  The  depth  scale  is 
derived  from  the  UE6E  velocity  function.  The  continuous,  light  line 
is  from  the  gravity  model  shown  in  figure  8.  The  short,  heavy  lines 
mark  the  presence  of  coherent  reflectors.  Notice  the  agreement  of 
the  gravity  model  with  coherent  reflections  between  vibration  points 
20  to  60  and  70  to  150. 
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Fig.  20.  An  interpretation  of  the  first  1.5  seconds  of  two-way 
travel  time  for  the  E-3  profile  shown  in  figure  13.  The  depth  scale 
is  derived  from  the  UE6E  velocity  function.  The  continuous,  light 
line  is  from  the  gravity  model  shown  in  figure  8.  The  short,  heavy 
lines  mark  the  presence  of  coherent  reflectors.  In  the  interval 
between  vibration  points  130  and  160  the  agreement  of  the  gravity 
model  with  the  reflectors  is  not  high.  The  resolution  in  this  area 
is  not  as  high  in  the  three-dimensional  gravity  model  as  in  the 
two-dimensional  model  presented  in  Goforth,  et.  al.  (1979). 


Fig.  21.  An  interpretation  of  the  first  1.2  seconds  of  two-way 
travel  time  for  the  N-3  profile  shown  in  figure  14.  The  depth  scale 
is  dervied  from  the  UE6E  velocity  function.  The  continuous,  light 
line  is  from  the  gravity  model  shown  in  figure  8.  The  short,  heavy 
lines  mark  the  presence  of  coherent  reflectors.  In  figure  18  it  is 
evident  that  line  N-3  runs  parallel  to  and  on  the  margin  of  a  synform. 
This  geometry  causes  the  reflectors  to  migrate  out  of  the  plane  of 
the  profile,  a  difficulty  which  is  impossible  to  overcome  quanti¬ 
tatively.  In  spite  of  this  problem  the  section  shows  more  detail 
below  the  basin  than  either  of  the  east  to  west  profiles 
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time.  The  volcanic  units  dip  to  the  west  at  approximatly  45°.  This  is 
due  to  tectonic  rotation  by  the  Carpetbag  and  Yucca  Fault  systems. 

Dip  slip  motion  on  the  Carpetbag  Fault  probably  postdates  the 
deposition  of  the  volcanics  and  hence  is  less  than  about  11  million 
years  old.  The  early  opening  of  the  basin  at  Yucca  Flat  is  probably 
part  of  the  general  east-west  extension  of  the  Great  Basin  associated 
with  the  progressive  development  of  the  San  Andreas  Fault  (Atwater, 
1970).  The  strike  slip  motion  results  from  a  progressive  reorientation 
of  the  stress  field  due  to  major  plate  interactions.  The  Yucca  fault 
represents  an  eastward  migration  of  the  locus  of  faulting  as  strike 
slip  motion  and  constraints,  such  as  the  Climax  Stock,  have  become  more 
important  (Oldow,  personal  communication,  1981).  The  Yucca  fault 
generally  has  less  than  200  meters  of  displacement  on  it.  The  saddle 
region,  which  was  not  traversed  by  seismic  lines  in  the  eastern  basin, 
is  likely  to  be  a  remnant  of  a  previous  topographic  high  related  to 
northwest  to  southeast  trending  faults  which  predate  the  volcanics.  In 
the  saddle  the  Yucca-Carpetbag  Fault  system  breaks  up  into  a  very  dif¬ 
fuse  zone  of  en  echelon  faults  known  as  the  Topgallant  system.  The 
geometrical  features  observed  in  the  Yucca  Flat  geophysical  model  now 
have  a  structural  interpretation  consistent  with  the  regional  tectonic 
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CONCLUSIONS 

The  application  of  a  systematic  gravity  inversion  procedure,  which 
utilizes  known  structural  and  density  information,  has  resulted  in  a 
model  of  the  Tertiary-Paleozoic  interface  at  Yucca  Flat,  Nevada.  This 
model  is  accurate  to  better  than  100  n  in  most  areas  and  is  capable  of 
predicting  basin  depths  to  a  similar  accuracy.  This  is  achieved 
without  recourse  to  complex  density  variations,  but  with  only  a  single 
density  contrast.  The  estimation  of  the  regional  anomaly  is  crucial  to 
the  interpretation,  and  an  effective  procedure  for  this  purpose  has 
been  demonstrated. 

The  seismic  reflection  method  is  useful  for  mapping  the 
Tertiary-Paleozoic  interface  in  areas  where  reflections  are  not 
obscured  by  high  velocity  layers,  such  as  welded  tuffs,  and  where 
three-dimensional  structural  effects  are  not  too  serious.  The  high 
attenuation  resulting  from  Quarternary  sediments  in  this  locality  means 
that  low  frequency  sources  and  seismometers  should  be  employed. 
VIBROSEIS®  techniques  using  low  frequency  sweeps  were  particularly 
suitable,  however,  resolution  was  correspondingly  low.  In  addition, 
automatic  velocity  determinations  are  suspect,  at  least  in  the  deeper 
areas,  because  of  structural  complications.  Processing  with  migration 
before  stacking  and  the  use  of  interval  velocity  determinations  in 
wells  helped  alleviate  this  difficulty.  Faulting  is  usually  well 
defined  on  the  seismic  reflection  profiles. 


The  combined  interpretation  of  geologic,  gravity  and  seismic  data 
depicts  the  general  structural  configuration  of  the  basin  very  well. 
Major  normal  faults  on  the  west  rotate  the  basin  down  about  a  broad 
hinge  zone  at  the  eastern  margin.  This  structural  development  occurred 
after  depostion  of  the  volcanics.  A  pre-existing  northwest  structural 
lineation  forms  a  saddle  which  divides  the  basin  into  north  and  south 
parts.  This  saddle  affects  later  structural  elements  such  as  the 
Carpetbag  and  Yucca  Fault  systems  in  the  central  saddle  region.  A  tec¬ 
tonic  pattern  of  increasing  strike  slip  motion  on  the  western  bounding 
faults  in  recent  time  can  explain  the  shift  of  activity  from  the 
Carpetbag  to  the  Yucca  Fault  systems. 

Body  wave  magnitude  variations,  when  properly  normalized  to 
constant  explosion  yield,  show  a  variation  which  correlates  with  the 
model  presented  here.  In  particular,  location  of  the  explosions  rela¬ 
tive  to  the  Carpetbag  fault  and  to  the  saddle  is  highly  significant. 

The  theoretical  response  of  this  structure  will  be  calculated  to 
further  verify  this  hypothesis  and  provide  a  predictive  capability. 
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PART  II 
INTRODUCTION 

Body  wave  magnitude  estimates  for  a  number  of  nuclear  explosions 
at  Yucca  Flat,  Nevada  have  been  analyzed  for  spatial  variation 
(Alevine,  personal  communication).  Each  explosion  was  corrected  to 
a  common  datum  by  the  use  of  well  established  scaling  laws  for  yield 
and  depth  of  burial  (Mueller  and  Murphy,  1971).  These  data,  when 
plotted  on  the  map  of  Yucca  Flat  in  figure  22,  show  a  progressive 
decline  of  about  0.5  magnitude  units  from  the  middle  of  the  valley 
toward  the  east  over  a  distance  of  four  kilometers.  This  variation 
complicates  the  analysis  and  calibration  of  seismic  yield  estimation 
and  source  discrimination  techniques. 

The  magnitude  variation  is  correlative  with  the  structure 
contours  of  the  top  of  the  Paleozoic  units  as  presented  in  the  FIB 
model  of  Yucca  Flat.  The  systematic  nature  of  the  variation  suggests 
a  deterministic  explanation  of  the  phenomena  in  terms  of  scattering 
of  the  outgoing  body  waves  by  the  local  structure. 

This  problem  is  not  soluble  in  any  simple  form  for  realistic 
geophysical  models  of  Yucca  Flat.  A  numerical  solution  might  be 
computed  by  any  of  several  techniques  such  as  finite  difference, 
finite  element  or  some  form  of  ray  tracing.  After  a  review  of  the 
literature  on  numerical  wave  propagation  techniques,  the  above 
alternatives  were  rejected  in  favor  of  the  Aki  and  Lamer  (1970) 


Fig.  22.  Contours  of  body  wave  magnitude  variation  (m^)  super¬ 
imposed  on  a  contour  map  of  the  Tertiary-Paleozoic  Interface  for 
Yucca  Flat,  Nevada.  The  m^  contours  were  modified  from  Alewine 
(personal  communication) . 
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collocation  method.  This  technique  promised  to  have  high  accuracy 
and  reasonable  computation  speed.  It  is  not  limited  to  plain  strain 
theory  like  the  finite  difference  calculation  or  high  frequency  like 
the  asymtotic  ray  theory,  also  the  solution  is  not  restricted  to  a 
limited  set  of  rays  or  modes. 
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DERIVATION  OF  THE  METHOD 


In  this  section  a  technique  introduced  into  seismology  by 
Aki  and  Lamer  (1970),  Lamer  (1970),  Bouchon  (1973)  and  Bouchon  and 
Aki  (1977),  will  be  developed  in  a  general  form.  The  method  assumes 
a  solution  in  each  homogeneous  layer  or  region  in  terms  of  a  series 
of  plane  waves.  The  boundary  conditions  between  the  regions  are 
matched  at  a  finite  number  of  points  on  the  boundaries.  In  numerical 
analysis  this  type  of  calculation  is  called  an  orthogonal  basis 
boundary  value  collocation  scheme.  The  solution  is  expanded  in  an 
orthogonal  basis  of  known  solutions  to  the  partial  differential 
equation  and  forced  to  satisfy  the  boundary  conditions  at  a  finite 
number  of  points.  The  interpolatory  characteristics  of  the  solution 
are  used  to  provide  an  approximate  satisfaction  of  the  boundary 
conditions  at  unsampled  boundary  points.  Due  to  certain  requirements 
of  the  reciprocal  method  that  will  be  used  to  apply  the  explosion 
source,  it  is  more  convenient  to  develop  the  method  in  terms  of 
potentials.  From  Lame's  theorem  (Aki  and  Richards,  1980,  p  68),  the 
equations  of  motion  for  a  homogeneous  elastic  medium  will  yield 
solutions  in  terms  of  a  scalar  and  a  vector  potential. 

la)  U  * 

lb)  v-<f  *  0  , 
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*0  jif  •  «’  Va 

Id)  f  -  V‘  f  . 

These  equations  include  separate  wave  equations  for  compressional , 
lc,  and  shear,  Id,  waves.  In  general  the  potentials  are  ambiguously 
determined  and  hence  are  useful  only  when  certain  symmetries  are 
present  in  the  formulation. 

Plane  wave  solutions  with  harmonic  time  dependence, 

2)  f  --  A  exp  +  kx  +  i  te))> 

will  be  assumed.  In  the  case  of  a  wave  Incident  on  a  horizontal 
interface  all  scattered  wave  vectors  will  be  contained  in  the  same 
vertical  plane  as  the  incident  wave.  An  alternative  statement  is 
that  k  and  t)  are  constant  for  all  waves  satisfying  the  boundary 
conditions.  Thus  the  coordinates  may  be  rotated  to  eliminate  one 
wavenumber  component  and  the  potentials  reduced  to  three  scalar 
potentials,  with  horizontally  polarized  shear  (SH)  motion  uncoupled 
from  compressional  (P)  and  vertically  polarized  shear  (SV)  motion 
(Aki  and  Richards,  1980,  p  215).  In  the  case  of  an  interface  with 
topography  variable  in  the  x-direction  only  and  incident  waves  making 
the  angle  SI  with  that  direction,  the  y  wavenumber  component,  t)  ,  will 
remain  constant,  but  the  x  wavenumber  component,  k  ,  will  couple  over 
the  entire  wavenumber  spectral  range.  The  constancy  of  t)  permits  the 
specification  of  a  vertical  plane  for  each  value  of  k  which  will  con¬ 
tain  the  scattered  waves.  The  P-SV  and  SH  motions  will  uncouple  for 
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each  k  and  only  three  scalar  potentials  are  required  (Larner,  1970, 
p  25-44).  The  new  propagation  coordinates  are: 


X  * 

30 

IXx'l  * 

I  kx  +  yy  1  , 

30 

K  ■ 

X  cm  Cn) , 
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/x’\ 

B  i  =  ( v ’)  = 
\2  j 
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The  displacements  are  written  in  terras  of  the  required  potentials  as: 


4a) 

u, 

<*’>  y\*) 

_  _ 
&<• 

d  2 

4S) 

Uj 

0*',  y\  ?•) 

A3L 

2)71  , 

4c) 

U, 

&  4 

-  -■?■■■ g-  -r 

d  a 

JliL 

dx' 

83 


In  each  homogeneous  region  Che  general  solution  will  be  expanded  in 
an  infinite  number  of  plane  wave  contributions: 

5^)  fi  Cx.v.i.k)  3  [a.W  ckP  (-,oO  ♦  Az  (y)  exp  0'v)i)J 

5b)  tx,y,i,k)  3  [A3  (w)  «*P  (-i-rz)  "A  <  (W)  exp  exp 

5c)XCx,y^,k)  =  | /\  5  (k)  exp(-i>i)  -  A„  (k)exp  (»>*)]  <t*p(->(i<x*9y)) . 

the  term,  exp(-iut),  is  suppressed  in  all  the  following  equations. 

The  boundary  conditions  to  be  satisfied  are  continuity  of 
displacement  and  traction  normal  to  the  boundary  at  each  internal 
interface.  The  normal  traction  is  set  to  zero  on  the  free  surface. 
The  general  solutions  for  the  displacement  and  stress  components  are: 


<~) 

U;  Cx.y.l)  * 

A*  exp  (,-ikx)  dK, 

Cob) 

Oij  C.*>y>*)  -- 

Ax  exp  (-',kx)  dk, 

' 

1,3, 

*  1 t  •••)(*  » 

The  normal  tractions  at  an  interface  are  obtained  from  the  inner 
product  of  the  unit  normal  vector  and  the  stress  tensor  or  the 
operator , 

1)  p,i  1  fij,  «j  (6)  .  i.j  *  '.*.»> 

The  unit  normal  is  a  function  of  the  interface  topography,  |  (x).  . 
The  operators  e,  f  and  F  are  exponential  functions  of  y,  z,  k  and  t)  , 
as  well  as  the  elastic  constants  and  density.  They  are  wavenumber 
domain  spacial  derivative  operators  as  required  by  4.  The  Internal 


boundary  conditions  require  that  these  quantities  be  evaluated  on 
each  side  of  the  Interface.  In  matrix  form  this  is  expressed; 


84 


exp  (-ikx)  dk  . 


The  integration  over  k  must  be  performed  numerically,  so  the  Fourier 
transform  is  replaced  by  a  Fourier  series  of  finite  length  (rectangle 
rule  integration  or  discrete  Fourier  transform).  This  results  in  a 
periodic  repetition  of  the  interface  topography  with  fundamental 
wavelength,  L.  The  collocation  technique  evaluates  the  displacements 
and  tractions  at  a  finite  number  of  points  in  the  interval  [o,l]. 

The  space  and  spectral  discretizations  are: 


O' 

*r> 

Ax  h  , 

n  =  o,i)  « • 

•i  10-/, 

krr> 
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♦  fa-0 
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At  this  point  the  notation  required  becomes  rather  complex.  In 
addition  to  the  tensor  indices  we  now  have  indices  associated  with 
the  series  summation  and  collocation  points.  The  differing  layers 
of  the  model  must  also  be  distinguished.  Figure  23  is  introduced  to 
clarify  some  of  the  additional  notation.  The  superscripts,  i,  will 
serve  to  identify  interfaces  starting  with  the  free  surface,  1*1. 

The  superscripts,  1,  refer  to  the  layers,  where  layer  1  is  just  below 


Fig.  23.  An  introduction  to  the  indicial  notation  used  to 
identify  the  homogeneous  regions  and  interfaces  used  in  the  Aki- 
Larner  technique  derivation. 


the  1th  interface.  There  are  P  layers  and  interfaces.  For  the 
discrete  space  specified  by  9,  equation  8  may  be  written. 
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id)  U  ‘  =  G‘  K . 

In  10,  U*  is  the  6*N  vector  of  displacements  and  tractions  evaluated 
at  the  discrete  interface  points,.  G*  is  the  6  •  Mx  6  •  H  matrix 
which  transforms  the  discrete  k,  6  •  M  vector  of  potential 
coefficients,  A*,  in  layer  i.  In  the  simplest  application  we  will 
set  M-N  for  a  square  linear  system  of  equations.  A  similar  equation 
exists  in  layer  i  +  1  at  the  same  interface,  which  would  be  the 
(i  +  l)th  interface.  If  the  and  Ui+*  vectors  evaluated  at 
^ (x)  are  set  equal  the  result  may  be  solved  for  A*  in  terms  of 


,0  A'  •  [S']"  &'*’  Au\ 

The  product  matrix  is  the  propagator  matrix  for  this  problem.  The 
product  of  each  such  propagator  at  all  Interfaces, 

«)  X  =  &’  W)  TT  [&■'  _ 

i-7 

results  in  a  relationship  between  the  displacement  and  stress  at  the 
free  surface  and  the  wave  field  in  the  half-space,  layer  P,  given  by 

>i)  U’  -  J  Ap 

An  incident  body  wave  is  specified  in  the  half  space,  then  the  free 
surface  boundary  conditions  of  vanishing  stress  are  used  to  solve  for 

p 

the  reflected  waves  in  the  halfspace  thus  completely  specifying  A  . 
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The  A*  vectors  in  each  other  layer  are  found  by  application  of  the 
propagator  matrices  in  equation  11. 

Larner  (1970)  introduced  a  variation  of  the  collocation  method 
which  is  of  considerable  importance.  The  partitions  of  U  and  the 
columns  of  the  partitions  of  G  are  Fourier  transformed  over  N  space 
samples.  The  boundary  conditions  are  then  specified  in  the  wave- 
number  domain.  The  resulting  Fourier  series  is  then  truncated  to  MSN 
points  and  the  equations  solved  as  before,  except  for  an  inverse 
transform  to  convert  the  displacements  and  tractions  back  to  the 
space  domain.  For  M=N  this  procedure  is  identical  to  the  space 
domain  collocation,  otherwise  it  is  a  least  squares  solution  of  the 
collocation  equation. 

The  wavenumber  spectrum  for  this  problem  possesses  poles 
corresponding  to  surface  wave  modes.  Integration  along  the  real 
wavenumber  axis  will  result  in  numerical  instability.  The  inclusion 
of  viscoelastic  attenuation  will  remove  the  poles  from  the  real 
k  -  axis  and  permit  integration  along  contours  parallel  to  that  axis. 
The  attenuation  may  be -specif  led  as  an  average  temporal  Q  for  the 
model  with  a  complex  frequency  and  wavenumber  or  as  an  average 
spatial  Q  with  complex  velocity  and  wavenumber.  In  order  that  the 
proper  phase  relationships  exist  across  the  incident  wavefront  the 
imaginary  part  of  the  horizontal  wavenumber,  in  all  layers,  must 
equal  the  Imaginary  part  of  the  horizontal  wavenumber  of  the  incident 
wave  (Aki  and  Larner,  1970).  The  imaginary  part  of  the  frequency  is 


determined  by, 

K)  ZLr~  C fS)  =  •  2tr  If  |  /  (2  Q)  , 

for  a  temporal  Q. 

The  numerical  accuracy  of  the  solution  technique  is  very  high. 
If  the  interface  topography  is  adequately  sampled,  so  that  no 
aliasing  occurs  (i.e.,  the  Fourier  expansions  are  convergent),  then 
the  solutions  are  computed  with  accuracy  equivalent  to  the  rounding 
error  of  the  computer.  The  simultaneous  equations  in  11  and  13  are 
usually  well  conditioned,  so  that  little  loss  of  significance 
occurs  in  their  solution.  A  different  form  of  error  results  from 
the  so  called  Rayleigh  ansatz  (Aki  and  Lamer,  1970).  For  interface 
slopes  greater  than  about  60°,  the  solution  should  contain  multiply 
reflected  waves  which  are  not  provided  for  in  the  specification.  In 
practice  this  error  has  not  been  very  important. 
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SOURCE-RECEIVER  RECIPROCITY 

The  Aki-Lamer  formulation  may  be  used  to  calculate  the 

response  of  an  irregularly  layered  model  to  an  incident  plane  body 

wave.  The  theory  must  be  extended  to  include  point  explosive 

sources.  Bouchon  and  Aki  (1977)  demonstrate  how  a  discrete  wave- 

number  representation  can  be  used  to  calculate  the  near  field 

response  of  complex  sources  in  irregularly  layered  structures. 

Because  of  our  interest  in  teleseismic  measurements  and  the  far  field 

response,  a  more  efficient  procedure  results  from  the  use  of  the 

seismic  reciprocity  theorem  (Bouchon,  1976)  and  the  plane  wave 

response.  The  reciprocity  theorem,  after  White  (1965),  states: 

If,  in  a  bounded,  inhomogeneous,  anisotropic  elastic  medium, 
a  transient  force  f(t)  applied  in  some  particular  direction 
at  some  point  P  creates  at  a  second  point  Q  a  transient 
displacement  whose  component  in  some  direction  is  u(t), 
then  the  application  of  the  same  force  f(t)  at  point  Q  in 
direction  will  cause  a  displacement  at  point  P  whose 
component  in  the  direction  is  u(t). 

The  response  of  a  point  dilatational  source,  located  at  P  in 

layer  j ,  at  a  point  Q,  in  the  half  space  layer  n,  can  be  derived  from 

the  consideration  of  a  force  F  at  Q.  The  force,  F,  is  directed  along 

the  ray  connecting  P  and  Q  and  results  in  a  vector  displacement  U(P) 

at  P.  By  the  reciprocity  theorem  the  same  force  at  P,  in  the 

displacement  direction,  will  produce  a  displacement  U(Q)  in  the  ray 

direction  (P-wave  motion).  The  dilatational  source  is  modeled  by 

three  orthogonal  force  dipoles.  To  obtain  the  displacement  at  Q,  in 


A  THEORETICAL  RESPONSE  PROFILE 


The  procedure  outlined  in  the  previous  sections  has  been  coded 
in  FORTRAN  IV.  The  program  has  been  tested  on  a  number  of  simple 
problems  including  flat  layered  and  gently  perturbed  structures. 
Working  versions  exist  for  both  the  Control  Data  Corporation  6600 
and  the  Digital  Equipment  Corporation  VAX  11-780. 

The  scattering  hypothesis  can  be  tested  on  a  profile  taken 
from  the  Yucca  Flat  model.  The  first  calculation  has  been  done  on 
an  east-west  profile  at  a  latitude  of  37°  02'  N.  The  basin  has 
very  little  variation  parallel  to  strike  and  good  structural  control 
is  provided  by  geophysical  surveys  in  that  area.  In  figure  24,  the 
observed  variation  south  of  37°  05'  is  plotted  along  with  a 
profile  taken  from  Goforth,  et.  al.  (1979).  The  calculated,  one 
Hertz  P-wave  amplitude  variation  for  sources  distributed  along  the 
profile  in  figure  24  should  agree  with  this  profile  if  the  hypothesis 
is  correct. 

The  geologic  structure  of  the  profile  can  be  approximated  by 
four  layers.  The  Paleozoic  limestone  can  be  modeled  as  a  half  space 
which  is  faulted  at  a  low  angle  (~ 50°)  on  the  west  by  the  Carpetbag 
fault  system  so  that  an  asymmetrical  basin  is  formed.  The  basin  is 
filled  with  Tertiary  volcanics  and  Quaternary  alluvium.  The  water 
table  at  a  depth  of  about  500  meters  divides  the  volcanics  into 
saturated  and  dry  units.  The  alluvium  is  entirely  free  of  water. 


Fig.  24.  A  profile  across  Yucca  Flat,  Nevada  at  a  latitude  of 
37°  02'  N  showing  the  observed  variation.  The  geological  struc¬ 
ture  was  taken  from  Goforth,  et.  al  (1979). 
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Average  velocities  and  densities  for  these  units  were  derived  from 
well  logs  as  contained  in  Table  III.  The  profile  is  shown  at  the 
bottom  of  figure  25. 

The  response  due  to  sources  at  the  four  locations  shown  in 
figure  25,  for  P-wave  angles  of  emergence  between  -40  and  40  degrees 
(west  and  east)  and  propagation  vectors  within  the  model  plane,  were 
computed.  An  average  Q  of  50  was  assumed  at  a  frequency  of  one 
Hertz  and  a  source  depth  of  0.8  kilometer  was  utilized.  This 
calculation  resulted  in  a  polar  radiation  pattern  of  amplitude  vs. 
emergence  angle  for  each  shot  point  as  shown  at  the  top  of  figure  25. 

In  order  to  model  the  averaging  process  inherent  in  the 
magnitude  calculation  the  median  amplitude  of  each  shot  will  be 
compared.  The  median  is  plotted  as  a  semi-circular  arc  on  each 
radiation  pattern.  The  ratio  of  median  amplitudes  between  shots  on 
the  west  (near  the  Carpetbag  fault)  and  those  on  the  east  is  2.4. 

This  corresponds  to  an  m^  variation  of  0.4  magnitude  units,  in  good 
ag'-eement  with  the  profile  in  figure  24. 

The  results  of  figure  25  are  limited  in  several  Important  ways. 
No  out  of  plane  rays  are  considered  and  only  a  single  narrow  frequency 
band  was  studied.  Future  calculations  will  include  the  out  of  plane 
rays  and  a  time  domain  synthesis  over  a  broader  (one  or  two  octave) 
frequency  band.  The  radiation  patterns  display  a  rather  complicated 
variation  which  may  be  observable  in  teleseismic  observations.  This 
phenomena  can  only  be  studied  by  examination  of  a  number  of  specific 
shots  and  receivers,  with  spectral  and  time  domain  modeling  of  the 


observed  P-waves. 


Fig.  25.  Polar  radiaCion  diagrams  of  teleseismic  P-wave  amplitude 
as  a  function  of  emergence  angle  for  shots  distributed  across  Yucca 
Flat,  Nevada. 
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CONCLUSIONS 

The  results  of  this  investigation  indicate  that  it  is  possible 
to  model  the  effects  of  near-source  scattering  on  teleseismic  signals 
from  Yucca  Flat,  Nevada.  This  is  made  possible  by  the  existence  of  a 
reasonably  accurate  geophysical  model  of  the  geologic  structure 
developed  from  gravity,  borehole  and  seismic  exploration.  The  Aki- 
Lamer  technique  is  a  computationally  effective  means  of  performing 
response  calculations  for  quite  complicated  models,  such  as  the  one 
exhibited  here.  Future  efforts  will  concentrate  on  time  domain  wave¬ 
form  synthesis  and  detailed  source  receiver  pair  models.  It  is  also 
possible  to  model  other  profiles  across  the  valley.  Similar  studies 
on  other  test  sites  can  guide  the  improvement  of  earthquake-explosion 
discrimination  and  yield  estimation. 
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LAJITAS  SEISMIC  STATION 


EUGENE  HERRIN 
GEOPHYSICAL  LABORATORY 
SOUTHERN  METHODIST  UNIVERSITY 


LAJITAS  SEISMIC  STATION 
Site  Report 

Introduction 

A  number  of  unpublished  seismic  noise  studies  during  the 
last  20  Years,  mostly  by  E.  Herrin,  W.  Guyton,  J.  Waugh  and 
B.  Brooks,  have  shown  that  sites  in  the  limestone  desert  region 
just  west  of  Big  Bend  National  Park,  Brewster  County,  Texas, 
show  very  low  short-period  noise  levels.  Peak-to-peak  displace¬ 
ments  at  the  surface  in  the  absence  of  wind  or  other  obvious 

-9 

disturbances  have  amplitudes  less  than  10  meter  at  a  frequency 
of  1  Hz.  and  less  than  10  meter  at  10  Hz. 

The  area  studied  is  very  sparsly  settled  and  has  only  one 
paved  highway  which  is  lightly  travelled.  The  nearest  railroad 
is  40  miles  to  the  west  and  has  only  a  few  trains  per  week. 

The  nearest  heavily  travelled  highway  and  railroad,  and  the 
nearest  town  of  any  size  -  Alpine,  Texas  -  is  about  75  miles  to 
the  north  of  the  area.  We  have  obtained  a  lease  on  private  land 
in  this  area  and  have  established  the  Lajitas  Seismic  Station 
which  is  described  in  this  report. 


Location 


The  Lajitas  Station  is  located  about  10  miles  by  road 
northeast  of  Lajitas,  Texas,  a  small  village  on  the  Rio  Grande 
just  west  of  Big  Bend  National  Park.  The  site  location  is 
shown  as  a  circled  dot  just  north  of  Terlingua  Sinkhole  on 
Figure  1,  and  is  on  land  owned  by  Mr.  Glen  Pepper.  Pertinent 
information  on  the  location  is  as  follows: 

Coordinates:  29  29'  02" 

103  40'  01" 

Elevation:  3325  feet  above  mean  sea  level 

The  paved  road  (Texas  FM  170)  and  graded  road  to  Pepper's 
hacienda  about  h  mile  south-southeast  of  the  site  can  be  seen 
in  Figure  1.  These  roads  provide  all-weather  access  to  the 
site  and  to  Lajitas,  Texas,  to  the  west  and  Study  Butte,  Texas, 
to  the  east.  The  remoteness  of  the  site  is  important  in  that 
man-made  noise  is  minimal  but  leads  to  certain  loaistical 
problems.  For  example,  the  nearest  scheduled  airline  service 
is  to  be  found  at  Midland,  Texas,  250  miles  to  the  northeast 
and  El  Paso,  Texas,  300  miles  to  the  northwest. 
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Figure  1 
Portion  of  map 

Amarilla  Mountain  Quadrangle 
Texas  -  Brewster  County 

7.5  Minute  Sevies  . 

Scale  1:  24,000  Contour  Interval  40  feet 


Geology 


The  Lajitas  site  lies  on  the  top  of  the  Santa  Elena  for¬ 
mation,  a  massive  limestone  unit  of  Cretaceous  (Upper  Albian) 
age.  Cretaceous  formations  beneath  the  site  are  listed  under 
the  heading  "Big  Bend  National  Park  Area"  in  the  correlation 
table  shown  in  Figure  2.  Table  1  gives  thicknesses  and  brief 
descriptions  of  the  underlying  Cretaceous  formations  as  seen 
in  exposures  in  Santa  Elena  Canyon  (SEC)  about  12  miles  south¬ 
east  and  The  Solitario  (SOL)  about  8  miles  northwest  of  the  site. 
The  Santa  Elena  Canyon  data  are  from  Maxwell  et.  al.,  1967, 

Univ.  of  Texas  Pubic.  6711,  3uv.  of  Economic  Geology,  and  The 
Solitario  data  are  from  Herrin,  1958,  Geology  of  the  Solitario 
area,  Trans-Pecos  Texas,  Harvard  University. 

From  Table  1  it  can  be  seen  that  the  site  is  underlain 
by  more  than  3600  feet  of  lower  Cretaceous  sediments  consisting 
predominantly  of  massive  limestones  which  should  have  P-wave 
velocities  in  excess  of  5  km/sec. 

The  uppermost  800  ft  consists  of  the  massive,  dense  Santa 
Elena  formation.  Beneath  the  basal  Cretaceous  conglomerate 
(Shut  Up  formation)  there  are  several  thousand  feet  of  folded 
Paleozoic  sediments,  mostly  siliceous  shales  with  some  chert  and 
limestone,  underlain  by  a  Cambrian  and  older  basement  complex. 

The  site  is  on  a  horst  bounded  on  the  east  by  a  normal 
fault  on  the  west  side  of  the  Long  Draw  graben  and  represented 
by  the  fault-line  scarp  seen  on  Figure  1  about  *s  mile  east  of 
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Table  1 

Thickness  (feet) 


Formation 

SEC 

SOL 

Description 

Santa  Elena 

740 

820 

Massive  hard  limestone 
with  bedded  chert 

Sue  Peaks 

265 

187 

Marly  limestone  -  Thin 
to  medium  bedded 

Del  Carmen 

465 

685 

Massive,  hard,  cherty 
rudistid  limestone 

Telephone  Canyon 

135 

190 

Nodular  limestone  and  marl- 
thin  to  medium  bedded 

Maxon 

0-10 

0 

Calcareous  sandstone 

Glen  Rose 

? 

1154 

Alternating  units  of  massive 
limestone  and  highly  fossili- 
ferous  marl 

Yucca 

7 

482 

Sandy  limestone  and  dolomite 
grading  into  massive  limestone 
at  top 

Shut  Up 

? 

100 

Chert  boulder  basal  conglomerate 
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Figure  2 

Correlation  Table  for  Cretaceous  Formations 
(from  Maxwell  et  al.,  1967,  Univ.  of  Texas 
Pubic.  6711,  Bur.  of  Economic  Geology) 


the  circled  dot.  The  horst  is  bounded  on  the  west  by  a  series 
of  normal  faults,  mostly  down-thrown  to  the  west,  culminating 
in  the  Well  Creek  graben  which  lies  just  north  of  California 
Hill.  It  is  possible  that  this  arrangement  of  normal  faults 
serves  to  shield  the  site  from  short  wave-length  surface  waves 
propagating  from  noise  sources  which  do  not  lie  on  the  horst. 
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Drilling  Operations 


Dick  Baker  Drilling  Company  of  Marfa,  Texas,  began  oper¬ 
ations  at  the  site  on  8  Sept.  1980.  Operations,  including  drilling, 
setting  casing  and  cementing  three  holes,  were  completed  on  19 
Sept.  1980.  The  three  holes  are  located  with  ten-foot  spacing  on 
a  north-northeast  bearing  line-of-centers .  The  north-most  hole  is 
labelled  number  1;  the  south-most  is  number  three.  Specifications 
for  the  holes  are  given  in  Table  2. 


Table  2 


(inches) 

(feet) 

Maximum 

Hole 

Diameter 

Total  Depth 

Inclination 

1 

9  7/8 

352 

2°  (at  350  ft) 

2 

7  :/8 

122 

1°  (at  120  ft) 

3 

9  7/8 

350 

1%°  (at  350  ft) 

Holes  number  1  and  3 

were  cased  and 

cemented  to  the  bottom 

with 

API  7"  OD  23#  casing. 

Hole  number  2 

was  cased  and  cemented 

to  the  bottom  with  API  5"  OD  N  80  casing. 

The  holes  were  air-drilled.  No  water  was  encountered  in 
drilling,  and  no  water  is  present  within  the  casings.  Nipples  were 
installed  so  that  1%  to  2  feet  of  casing  are  exposed  above  the  sur¬ 
face  with  male  threads  up  and  protected  with  suitable  caps. 
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Facilities 

A  rectangular  concrete  pad  was  poured  around  the  three  cased 
boreholes  extending  at  least  5  feet  horizontally  from  all  holes. 

The  site  has  been  graded  to  provide  a  parking  area  for  vehicles 
and  trailers,  and  the  road  to  the  site  from  Pepper's  main  road 
has  been  graded  and  filled.  The  site  is  easily  accessible  by 
passenger  car. 

An  adobe  hut  with  approximately  15  ft.  x  15  ft.  of  interior 
floor  space  has  been  constructed  adjacent  to  the  concrete  pad. 
Electrical  power  (110  and  220  volt,  single  phase)  was  brought  to 
the  site.  -Three  power  drops  were  placed  to  serve  the  hut  and  at 
least  two  trailers.  Telephone  lines  were  strung  (Big  Bend  Telephone 
Company)  to  provide  a  leased  data  line  to  Dallas  (SMU  campus)  and 
a  standard  telephone. 

A  water  line  was  laid  from  Pepper's  main  water  supply  to 
provide  water  for  use  in  trailers  and  for  fire  protection.  A 
drain  system  and  septic  tank  was  installed  so  that  a  camper- 
trailer  could  be  placed  on  the  site. 


TWJgPBEIJBP 
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INSTRUMENTATION 

Initial  tests  were  begun  in  December,  1981,  using  a 
Teledyne-Geotech  23900  vertical  seismometer  in  the  Number  2 
Lajitas  borehole  with  a  Teledyne-Geotech  42.21  amplifier 
and  resistive  damping.  This  short-period  system  was  designed 
to  give  approximately  a  flat-velocity  response  to  ground 
motion  between  1  Kz.  and  10  Hz.  A  series  of  tests  between 
January  and  June,  1981,  showed  that  the  resolution  of  this 
system  was  limited  by  system  noise  for  seismically  quiet 
periods  at  frequencies  above  about  5  Hz.  On  24  June,  1981, 
a  Geotech  "brick”  amplifier-prototype  GS  43310-  using 
feedback  damping  was  installed  at  the  Lajitas  site,  replacing 
the  42.21  amplifier.  With  this  new  configuration  it  was 
possible  to  resolve  "quiet-period"  ambient  background  noise 
to  frequencies  as  high  as  20  Hz.  The  response  of  the  system 
was  approximately  flat  to  ground  velocity  from  2  Hz.  to  17  Kz. 
with  the  response  down  6  dB  at  1  Hz.  and  3  dB  at  20  Hz. 

The  23900-  "brick"  system  remained  in  use  at  the  Lajitas 
site  until  Fall,  1982.  A  typical  "quiet-period"  displacement 
background  spectrum  for  Lajitas  is  shown  in  Figure  3  in  terms 
of  decibels-seismic  vs.  log  frequency.  It  can  be  seen  that 
the  spectral  level  is  limited  by  system  noise  rather  than 
ground  displacement  for  frequencies  above  about  20  Hz.  and 
displacement  power  below  about  -225  dBs .  This  system  noise 
level  was  predicted  by  O.D.  Starkey  (personal  communication) 
of  Teledyne-Geotech  based  on  the  design  specifications  of 
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the  seismometer  and  the  amplifier. 

In  June  of  1981  it  became  clear  that  a  more  remote  site 
was  needed  near  the  Lajitas  station.  Activity  going  on 
around  the  Lajitas  site  was  substantially  increasing  the 
ambient  seismic  background  noise.  Also  it  was  planned  that 
AFTAC  would  take  over  Lajitas  holes  numbers  1  and  2  for 
inclusion  in  their  SOUS  (Southern  U.S.)  tripartite  network. 
On  June  7,  1981,  a  remote  site  was  established  in  an 
abandoned  mine  on  Tres  Cuevas  mountain  about  3.6  miles 
(5.8  km)  WSVJ  of  the  Lajitas  site.  At  the  Tres  Cuevas  site 
we  used  a  Teledyne-Geotech  18300  short-period  vertical  seis¬ 
mometer  with  a  "brick"  amplifier.  The  frequency  response 
and  system  noise  levels  were  essentially  the  same  as  for 
the  Lajitas  bore-hole  system.  The  Tres  Cuevas  site  is 
remote  from  any  travelled  road^access  is  by  jeep-road  with 
locked  gates.  The  station  is  powered  by  solar  panels,  and 
the  data  is  transmitted  by  VHF  telemetry  to  the  Lajitas 
site.  "Quiet-period"  background  noise  spectra  at  the  two 
sites  are  essentially  the  same;  however,  the  mine  site  noise 
is  not  affected  by  activity  at  the  Lajitas  station. 

In  July  of  1981  Teledyne-Geotech  carried  out  sone 
preliminary  array  studies  at  the  Lajitas  site  and  in  July 
of  1932  Sandia  National  Laboratory  began  extensive  high- 
frequency  array  experiments  at  the  site.  These  experiments 
are  still  going  on. 

In  the  summer  and  fall  of  1982  AFTAC  completed  work  on 
their  SOUS  stations  installing  23900  seismometers  near 
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Marathon,  Texas;  Shatter,  Texas;  and  in  hole  number  2  at 
the  Lajitas  site.  A  KS  36000  system  was  installed  by  AFTAC 
in  Lajitas  hole  number  1.  Under  an  agreement  involving  the 
U.S.  Air  Force  (AFTAC),  DARPA  and  SMU  the  data  from  the 
SOUS  stations  were  made  available  to  us  at  the  Lajitas  site. 
We  are  currently  sending  short-period  seismic  data  by 
FM-telemetry  from  the  Marathon,  Shatter,  Lajitas  and  Tres 
Cuevas  (mine)  sites  back  to  our  laboratory  in  Dallas.  The 
locations  of  these  stations  in  the  Big  Bend  region  of  Texas 
are  shown  in  Figure  4 .  The  coordinates  of  the  stations 
are  given  below: 


Lajitas 

29 

20* 

02" 

n; 

103 

40' 

01 

Tres  CudVas 

29’ 

13’ 

59" 

h; 

103° 

43' 

05 

(mine) 

Shafter 

29° 

55’ 

27" 

n; 

104° 

22' 

17 

Marathon 

30° 

H* 

CO 

22" 

n; 

103’ 

15' 

17 

Figure  4 ♦  Map  of  Big  Bend  area  of  Texas  ^ 
showing  Tres  Cuevas,  Villa  de  la  Mina,  Shafter 
and  Marathon. 
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DATA  ANALYSIS 

Using  FM  telemetry  over  leased  telephone  lines,  data 
has  been  or  is  being  transmitted  from  the  SOUS  sites  (plus 
Tres  Cuevas)  -  see  Figure  4  -  as  well  as  from  KS  36000 
seismometers  at  Albuquerque  (ANMO-SRO)  and  McKinney,  Texas, 
to  our  laboratory  at  SMU  -  see  Figure  5 •  Recently  we  have 
elected  to  substitute  data  from  the  KS  36000  located  on 
campus  for  the  McKinney  data  and  so  eliminate  one  of 
the  FM  telemetry  links.  At  the  price  of  increased  ambient 
background  noise  we  have,  by  this  change,  provided  a  source 
of  high  dynamic  range  data  (not  limited  by  the  FM  channel) 
which  can  give  us  unclipped  digital  data  for  large  events 
which  will  cause  all  the  other  channels  to  be  clipped. 

Figure  6  shows  the  data  flow  from  the  field  sites 
to  Dallas  with  the  exception  of  the  SMU  campus  site  being 
substituted  for  the  McKinney,  Texas,  site. 

Data  from  the  FM  lines  are  routed  :.nto  analog-to- 
digital  conversion  system  included  in  a  DEC  11/23  dedicated 
computer  as  shown  in  Figure  7.  Detection  is  performed 
using  the  Walsh  detection  algorithm  (Goforth  and  Kerrin, 

Bull.  Seismol.  Soc.  Amer.,  vol.  71,  no.  4,  pp.  1351-1360, 
1981) .  A  detection  bulletin  is  printed  and  the  wave-forms 
containing  the  detected  events,  sampled  at  40  sps,  are 
transmitted  to  the  VAX  11/750  computer  for  off-line  analysis. 

Since  fall  of  1982  we  have  been  looking  at  events 
detected  at  the  Lajitas  site,  particularly  for  near-regional 


Figure  7  Digital  acquisition  and  detection  system  at  SMU  Geophysical  Laboratory. 
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to  regional  distance.  Analysis  of  teleseismic  events  will 
begin  in  the  next  contract  period.  Data  for  20  regional 
events  are  shown  in  Table  2.  Distances  were  estimated  from 
single  station  data  (Tres  Cuevas  site)  based  on  Pn  and  Lg 
times  (Pn  velocity  8.0  km/sec,  Lg  velocity  3.5  km/sec., 

Pn  intercept  time  8  sec) .  The  data  was  read  from  analog 
recordings,  and  corrections  were  made  for  instrument 
response  based  on  visual  estimation  of  the  predominant 
period.  Magnitudes  were  calculated  using  the  following 
expressions : 

A 

For  Pn:  Mb  =  -17.0  +  log  IT  +  7  log  D, 

A  in  millimicrons  (o-p)  and  distance  (D)  in  km. 

Based  on  Evernden's  results  (Bull.  Seismol.  Soc.  Amer., 
vol.  57,  n.  4,  p.  611,  Eqs.  7  and  8)  for  propagation 
of  Pn  from  Nevada  Test  Site  eastward  through  Texas. 

For  Lg;  Mb  =  -1.73  +  2.56  log  D  +  log  A 
with  200  km  <_  D  >_  800  km. 

A  is  maximum  sustained  amplitude  in  microns  (o-p) . 

Based  on  Mitchell  and  Nuttli ' s  results  for  Lg  in  Northern 
Iran  (Semi-Annual  Tech.  Report,  1  October  1977  -  31  March 
1978,  AFOSR  -  DARPA,  St.  Louis  University' . 

As  can  be  seen  in  Table  2,  the  estimated  magnitudes 
based  on  Lg  and  Pn  show  surprising  little  scatter  -  the 
average  value  was  taken  as  representative  of  the  event. 

The  signal  to  noise  ratio  was  estimated  for  the  Pn  and  Lg 
signals  associated  with  each  event.  This  ratio  was  normalized 


♦Slightly  clipped 


LAJITAS  -  REGIONAL  SIGNALS  -  TRES  CUEVA6  SITE 
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for  an  event  of  Mb  2.0.  The  S/N  in  dB  were  plotted  against 
log  distance  to  show  two  roughly  linear  trends ,  for  Pn  and 
Lg.  These  trends,  all  for  Mb  =  2,  were  as  follows: 

Lg :  35  dB  at  200  km  down  to  about  5  dB  at  600  km 

Pn:  28  dB  at  200  km  down  to  -  5  dB  at  600  km 

These  results  give  a  basis  for  determining  the  detection 
capability  of  a  single  seismometer  at  the  Latjias  site 
for  Pn  and  Lg  arrivals  at  regional  distances.  Inferences 
as  to  detection  capabilities  for  magnitudes  other  than  2.0 
can  be  made  using  simple  scaling  but  only  for  the  distance 
range  200  to  600  km. 

Examples  of  two  teleseismic  events  from  the  USSR  are 
shown  in  Figures  8  and  9.  It  should  be  noted  that  the 
local  (Lajitas)  magnitude  is  about  0.3  units  smaller  than 
the  NEIS  mean  magnitude.  Whether  this  difference  is  real 
or  a  result  of  network  bias~*in  the  UET5~ data  Ts“  not  known 
at  this  time.  Clearly  one  of  the  more  important  objectives 
of  our  on-going  research  program  is  to  determine  if  an 
exceptionally  quiet  site  like  Lajitas  also  produces  abnormally 
low  teleseismic  signal  amplitudes .  An  answer  to  this 
question  awaits  the  accumulation  of  data  using  the  automated 
seismic  analysis  system  now  in  development. 


Fiqure  8.  SEISMIC  TRACES  OF  THE  USSR  NUCLEAR  EXPLOSION  ON  MAY  25,  1981  AT  SEISMOMETERS  2  AND  H.  THE  ORIGIN  TIME  IS 
— '  04:59:57.2  UTC  AND  THE  LOCATION  IS  SOUTH  OF  NOVAYA  ZEMLYA  (68.182N.053.689EI  AT  A  DISTANCE  OF  76°  THE 

BODY  WAVE  MAGNITUDE  IS  5.5. 
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ABSTRACT 


The  resolution  of  the  seismic  sensors  discussed  in 
this  paper  is  limited  either  by  the  self-noise  of  the  system 
or  the  ambient  ground  noise  at  the  observing  site.  The 
resolution  limits  of  a  particular  system  are  frequency 
dependent.  A  major  objective  of  the  designer  is  to  insure 
that  seismic  sensors  are  available  which  will  be  limited 
in  resolution  only  by  the  ambient  background  noise  at  the 
quietest  sites  over  the  frequency  band  of  interest  in 
treaty  verification  research.  This  paper  presents  a 
review  of  the  data  currently  available  on  the  limiting 
resolution  of  the  most  advanced  instruments  used  in  this 
research  program. 


INTRODUCTION 


The  seismometer-amplifier  combinations  discussed  in 
this  paper  are  the  "state-of  -the-art"  systems  now  in  use  or 
in  development  for  use  in  treaty  verification  research.  Only 
systems  which  can  be  placed  in  bore-holes  with  casings  of 
7  inches  or  less  in  diameter  fall  into  this  category,  because 
the  bore-hole  environment  has  been  found  to  provide  the 
optimum  stability  and  lowest  ambient  background  noise  at 
any  given  site.  All  of  these  instruments  are  manufactured 
by  Teledyne-Geo tech  (TG)  in  Garland,  Texas,  thus  the 
Teledyne  model  numbers  will  be  used  throughout  this  paper  in 
the  discussions  of  the  properties  of  the  instruments.  Systems 
currently  in  use  are  the  TG  36000-01  in  the  Seismic  Research 
Observatories  (SRO) ,  the  TG  36000-04  and  TG  750  in  the 
National  Seismic  Systems  (NSS)  and  the  TG  36000-03  and  TG  20171 
(23900)  used  at  installations  operated  by  the  Air  Force  Technical 
Applications  Center  (AFTAC) .  The  TG  23900  is  a  version  of  the 
TG  20171  designed  for  operation  in  deep  bore-holes,  but  the 
two  are  essentially  the  same  instrument  in  so  far  as  the 
discussions  in  this  paper  are  concerned. 

Properties  of  The  Instruments 

The  instrumental  characteristics  of  most  importance  in 
treaty  verification  research  can  be  separated  into  four 
categories:  Effective  band-width,  dynamic  ranae ,  linearity 


at  about  -  220  dBs  as  the  frequency  approaches  10  Hz.  Is 
the  apparent  resolution  limit  at  about  -  220  dBs  being  set 
by  the  ambient  ground  motion  or  by  the  self-noise  of  the 
TG  36000-01? 

In  Figure  3  the  noise  spectral  densities  in  dBs  are 
shown  for  two  of  the  quieter  SRO's  (ANMO  in  New  Mexico  and 
BCAO  in  the  Central  African  Republic) .  The  BCAO  data  was 
taken  at  night  when  the  background  noise  at  that  station  is 
at  its  lowest  level.  Note  that  both  spectra  fall  sharply 
with  increasing  frequency  through  1  Hz.  The  "bumps”  in  the 
ANMO  spectrum  at  2  to  3  Hz  can  be  attributed  to  man-made 
noise  from  the  nearby  town  of  Albuquerque  and  from 
Kirtland  Air  Force  Base.  Beyond  4  Hz,  however,  the  two 
curves  approach  a  minimum  value  of  -210  to  -215  dBs. 

Data  such  as  that  shown  in  Figure  3  and  a  series  of  careful 
tests  carried  out  by  Sandia  National  Laboratory  (  H.  B.  Durham, 
personal  communications,  1982)  indicate  that  the  resolution 
of  the  TG  36000-01  is  set  by  self-noise  at  an  equivalent 
ambient  noise  level  of  -210  to  -220  dBs. 

This  conclusion  means  that  Peterson's  composite  noise 
model  (Figure  2)  does  not  represent  ambient  ground  noise 
at  frequencies  above  4  Hz  because  of  the  limited  resolution 
of  the  TG  36000.  In  order  to  measure  the  noise  levels  at 
sites  quieter  than  -220  dBs,  presumably  at  freauencies 
above  5  Hz,  we  must  use  an  instrument  with  lower  self-noise 
than  the  TG  36000. 
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measure  of  displacement  power  spectral  density  decibels 
relative  to  1  m  V  Hz.  which  we  call  decibels-seismic  (dBs) . 

The  effective  limit  of  resolution  of  seismometer- 

be. 

amplifier  combinations  canAset  by  (1)  the  maximum  gain  of 
the  system,  (2)  the  self-noise  of  the  system,  or  (3)  the 
ambient  background  displacement  at  the  operating  site. 

Clearly  for  best  performance  in  treaty  verificaiton  research, 
the  system  resolution  should  be  limited  by  the  ambient 
background  noise  over  the  entire  frequency  band  of  the 
system.  For  the  instruments  discussed  here,  the  maximum 
gain  is  great  enough  so  that  the  resolution  is  always 
set  by  either  the  self-noise  of  the  system  or  the 
ambient  background  noise. 

Seismic  Background  Noise 

Peterson  (1980)  reported  on  a  study  of  ambient  back¬ 
ground  noise  at  12  SRO  sites  (ANMO,  ANTO,  BCAO,  BOCO,  CHTO, 
GRFO,  GUMO,  MAIO,  NWAO,  SHIO,  SNZO  and  TATO) .  An  attempt 
was  made  to  determine  the  minimum  ambient  noise  at  each  site 
by  selecting  the  "quietest"  sample  from  the  digital  data 
available  for  the  study.  Using  these  "quiet  samples", 
displacement  power  spectral  densities  in  dBs  were  estimated 
for  each  site.  The  lowest  spectral  density  for  each  period, 
as  seen  at  any  of  the  12  sites,  was  then  selected  to  form 
a  composite  minimum  noise  model  for  these  sites  well 
distributed  over  the  earth.  Figure  1  shows  this  comoosite 
spectral  model  plotted  in  dBs  vs.  period.  The  individual 


points  shown  in  this  figure  indicate  the  spectral  estimates 
made  by  Fix  (1972)  using  data  recorded  at  the  Queen  Creek 
mine  aboet  50  km.  south  of  Pheonix,  Arizona.  Since  the  mid- 
1970’ s,  the  Queen  Creek  noise  spectrum  has  been  used  as  a 
low-noise  model  in  designing  seismometer-amplifier  combina¬ 
tions.  All  of  the  spectral  data  shown  in  this  paper  are 
for  vertical  seismic  sensors;  no  data  from  horizontal 
sensors  are  considered  here.  However,  the  conclusions 
reached  here  concerning  the  vertical  component  also  apply 
to  horizontal  components  of  the  TG  36000  and  TG  750  systems. 
The  TG  20171  has  only  the  vertical  sensor. 

One  conclusion  which  can  be  reached  from  Petersons' (1980) 
study  and  from  Figure  1,  is  that  for  periods  above  1  sec 
the  resolution  of  the  TG  36000-01  used  in  the  SRO's  is 
limited  by  the  ambient  seismic  background  noise  at  the  site 
and  not  by  the  self-noise  of  the  system.  This  conclusion 
also  applies  to  the  TG  36000-03  and  the  TG  36000-04  used 
in  AFTAC  and  NSS  installations  respectively,  again  for 
periods  above  1  sec . 

At  periods  less  than  1  sec  (frequencies  above  1  Hz) 
the  results  are  less  clear.  Peterson's  composite  low-noise 
spectrum  lies  below  the  Queen  Creek  values  as  might  be 
expected,  because  the  composite  spectrum  was  intended  to 
represent  a  lower  bound  on  ambient  background  noise.  This 
spectral  difference  can  be  better  seen  in  Figure  2  where  the 
data  from  Figure  1  have  been  replotted  as  dBs  vs.  frequency 
in  Hz.  Note  that  both  curves  in  Figure  2  tend  to  converoe 


and  resolution.  All  three  versions  of  the  TG  36000  are 
relatively  broad-band  instruments  designed  to  operate  over 
a  frequency  band  of  0.01  to  20  Hz.  The  TG  20171  and 
TG  750  are  short-period  systems  currently  being  used  with 
an  effective  band-width  of  about  0.5  to  40  Hz,  although 
in  some  installations  this  band -width,  for  various  reasons, 
is  truncated  by  low-pass  filters. 

All  of  the  systems  discussed  in  this  paper  have 
demonstrable  dynamic  range  in  excess  of  120  dB,  although 
because  of  the  configuration  and  over-all  frequence 
response  selected  by  the  user,  this  dynamic  range  may 
not  extend  over  the  entire  frequency  band. 

.  The  linearity  of  systems  with  such  a  large  dynamic 
range  is  very  difficult  to  measure  in  the  laboratory; 
however,  lower  bounds  can  be  set.  The  systems  discussed 
here  all  have  linearity  greater  than  70  dB  (about  one 
part  in  three  thousand) . 

The  last  category,  resolution,  is  the  primary 
subject  of  this  paper.  Hera  resolution  is  taken  to  be 
the  smallest  ground  displacement  that  can  be  resolved  by 
the  instrument.  Because  resolution  may  vary  with  frequency 
we  find  that  it  can  best  be  expressed  in  terms  of  dis¬ 
placement  power  spectral  density.  These  densities  are  in 
terms  of  true  ground  motion;  that  is,  system  outputs 
must  be  corrected  for  the  frequency-dependent  displacement 
response  of  the  seismometer-amplifier  combination  and 
any  subsequent  response  shaping  filters.  Following 
Peterson  (1980)  we  have  elected  to  use  as  a  standard 


A  New  Low  Noise  Model 


About  three  years  ago  the  author  began  a  series  of 
noise  measurements  at  a  very  remote  site  in  the  Dig  Bend 
region  of  Texas.  The  site,  located  near  the  small  village 
of  Lajitas  on  the  Rio  Grande  River,  is  at  29°  20'  north 

Q 

latitude  and  103  40'  west  longitude.  A  system  using  the 

TG  20171  vertical  seismometer  and  an  extremely  low-noise 
amplifier  (Teledyne-Geotech  43310)  was  used  to  measure 
the  ambient  background  noise  levels  at  the  Lajitas  site. 

Figure  4  shows  the  displacement  power  spectrum  for  Lajitas 
background  noise  along  with  the  Queen  Creek  spectrum.  Above 
about  2  Hz  the  Lajitas  spectrum  lies  about  10  dB  below 
the  Queen  Creek  values,  very  nearly  overlyincr  Peterson’s 
composite  low  noise  spectrum  out  to  about  4  Hz.  Beyond 
that  point  the  Lajitas  spectrum  continues  to  drop  and 
finally  flattens  at  about  -255  dBs  just  beyond  20  Hz, 
where  we  conclude  that  the  self-noise  of  the  TG  20171/43310 
system  begins  to  limit  the  resolution.  The  Lajitas  noise 
spectrum  is  now  being  used  as  a  low  noise  model  from  0.5  to 
20  Hz  in  the  design  of  new  seismometer-amplifier  systems. 

Conclusion  Regarding  Seismometer  Resolution 

Figure  5  is  a  diagram  prepared  by  O.  0.  Starkey 
(personal  communication,  1982)  showinq  the  resolution  limits 
of  the  systems  discussed  in  this  paper  as  set  by  the  calculated 
self-noise  of  the  seismometer-amplifier  combinations.  The 


Line  marked  "minimum  background  estimate"  is  the  new  low 
noise  model  based  on  the  Lajitas  noise  measurements. 

The  resolution  limit  for  the  SRO  and  AFTAC  broad-band 
instruments  (TG  36000-01  and  03)  is  close  to  -220  dBs 
which  roughly  corresponds  to  the  lowest  spectral  measure¬ 
ments  shown  in  Figures  2  and  3.  The  resolution  limit  for 
the  TG  20171/43310  used  to  measure  the  Lajitas  noise 
spectrum  is  about  -255  dBs  at  20  Hz.  This  limit  is  clearly 
indicated  in  the  Lajitas  noise  spectrum  shown  in  Figure  4. 
Resolution  limits  for  the  TG  36000-04  and  TG  750  used 
in  the  NSS  borehole  system  are  also  shown  in  dBs  as  a 
function  of 4 frequency .  Clearly  none  of  the  existing 
systems  are  capable  of  resolving  the  Lajitas  ambient 
background  noise  at  frequencies  above  20  Hz. 

A  new  system  being  developed  by  Teledyne-Geotech  for 
the  Defense  Advanced  Research  Projects  Agency  has  the  model 
number  TG  44000.  The  predicted  resolution  limit  for  this 
system,  as  shown  in  Figure  5,  is  about  -270  dBs  for  frequencies 
above  2  Hz.  If  the  TG  44000  system  now  beinq  tested  meets 
the  design-goals,  it  will  be  the  only  seismometer-amplifier 
system  which  can  resolve  the  ambient  ground  motion  at  the 
quietest  known  sites  over  the  frequency'  band  from  0.01  Hz  to 
beyond  20  Hz.  In  treaty  verification  research  we  are  con¬ 
stantly  searching  for  quieter  sites  at  which  to  make  seismic 
observations,  and  we  therefore  require  higher  resolution 
instruments  for  use  at  such  sites.  If  ambient  noise  measure¬ 
ments  are  made  using  instruments  with  resolution  limited  by 


their  own  self-noise,  we  can  not  determine  that  the  ambient 
noise  level  is  beneath  that  limit.  On  the  other  hand,  if 
we  do  not  have  the  measurements  to  indicate  that  a  site  is 
very  quiet  at  certain  frequencies,  we  may  lack  motivation 
for  developing  instruments  with  greater  resolving  power. 

In  the  'verification  research  program,  we  are  attempting  to 
develop  instruments  with  resolving  power  capable  of 
operating  at  the  quietest  sites,  with  a  resolution  limited 
only  by  the  ambient  seismic  noise  level. 
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FIGURE  CAPTIONS 


Figure  1. 


Figure  2. 


Figure  3. 


Figure  4. 


Figure  5. 


Composite  low-noise  spectrum  obtained  by  combining 
low  noise  points.  The  X's  are  for  Queen  Creek 
data  published  by  Fix  (1972)  (After  Peterson, 1980) 

Composite  low  noise  from  Figure  1  replotted  vs. 
frequency  (dots)  and  Queen  Creek  noise  spectrum 
(x's).  (After  Peterson,  1980,  and  Fix,  1972) 

Minimum  noise  spectra  for  the  SRO's  ANMO  (dots) 
and  BCAO  (x's)  (Modified  from  Peterson,  1980) 

Spectra  of  Lajitas  background  noise  (dots)  and 
Queen  Creek  noise  (x's). 

Limiting  displacement  resolution  of  a  number  of 
systems  based  on  the  calculated  self-noise 
levels  (Calculations  by  O.D.  Starkey,  1982). 
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Anomalous  Rayleigh  Waves  from  Nuclear 
Explosions  at  the  USSR  Shagan  River 
Test  Site 


by 

Tom  Goforth,  Bijan  Rafipour,  and  Eugene  Herrin 


ABSTRACT 


A  study  was  made  to  determine  if  signigicant  Dhase 
differences  exist  between  Rayleigh  waves  which  have 
traveled  almost  identical  paths  from  the  USSR  Shaaan  River 
test  site.  Surface  waves  from  five  nuclear  explosions 
recorded  at  six  SRO/ASRO  dicrital  stations  were  used  for 
the  analysis.  The  explosion  of  4  Aucrust  1979  was  selected 
as  a  reference,  and  at  each  site  the  spectral  Phases  of 
the  Rayleigh  waves  from  the  other  five  events  were  compared 
to  the  phase  of  the  reference.  A  technique  of  ohase- 
matched  filterino  (Herrin  and  Goforth,  1979)  was  used  to 
analyze  the  signals.  This  technique  reduces  the  effects 
of  multipathing  and  removes  phase  differences  due  to 
dispersion  along  slightly  different  travel  oaths. 

Each  epicenter  has  been  re-located  (North  and  Fitch, 
1980)  by  using  calibration  data  from  the  craterincr  shot  o* 

15  January  1965,  and  remaining  errors  in  location  and  oricrin 
time  are  considered  to  be  extremely  small.  Seismocrams  were 
analyzed  for  surface  waves  recorded  at  Matsushiro,  Japan 
(MAJO) ;  Shillong,  India  (SHIO) ;  Kabul,  Afghanistan  (KAAO) j 
Ankara,  Turkey  (ANTO) ?  Grafenburg,  W.  Germany  (GRFO) ;  and 
Alburquerque ,  N.M.  (AMMO) . 

Results  of  the  study  indicate  that  Rayleigh  waves 
from  some  Shagan  River  explosions  have  undergone  larcre 
ohase  shifts  relative  to  Rayleiah  waves  which  have  traveled 
almost  identical  paths  from  other  Shagan  River  explosions. 


In  some  cases ,  the  phase  shifts  can  be  interpreted  as  com¬ 
plete  phase  reversals  with  associated  time  delay.  In 
particular,  as  compared  to  the  explosion  of  4  August 
1979,  Rayleigh  waves  from  the  explosion  of  7  July  1979 
are  reversed  in  polarity  at  KAAO  and  ANTO  and  are  reversed 
in  polarity  and  delayed  at  GRFO,  SHIO,  and  MAJO. 
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